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Abstract 



Aerotaxis is the particular form of chemotaxis in which oxygen plays the role of 
both the attractant and the repellent. Aerotaxis occurs without methylation adapta- 
tion, and it leads to fast and complete aggregation toward the most favorable oxygen 
concentration. Biochemical pathways of aerotaxis remain largely elusive, however, 
aerotactic pattern formation is well documented. This allows mathematical modeling 
to test plausible hypotheses about the biochemical mechanisms. Our model demon- 
strates that assuming fast, non- methylation adaptation produces theoretical results 
that are consistent with experimental observations. We obtain analytical estimates 
for parameter values that are difficult to obtain experimentally. 
Chemotaxis in growth cones differs from gradient sensing in other animal cells, be- 
cause growth cones can change their attractive or repulsive response to the same 
chemical gradient based on their internal calcium or cAMP levels. We create two 
models describing different aspects of growth cone guidance. One model describes 
the internal switch that determines the direction of movement. However, this model 
allows chemotaxis under certain conditions only, so a second model is created to 
propose a mechanism that allows growth cone guidance in any environment. 
Endothelial cells go through extensive morphological changes when exposed to shear 
stress due to blood flow. These morphological changes are thought to be at least 
partially the result of mechanical signals, such as deformations, transmitted to the 
cell structures. Our model describes an endothelial cell as a network of viscoelastic 
Kelvin bodies with experimentally obtained parameters. Qualitative predictions of 
the model agree with experiments. 
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Chapter 1 



Introduction 



Biology has gone through an extraordinary change in the past century, partially due 
to increasingly advanced methods of being able to collect data, and partially be- 
cause of the sophistication in the quantitative analysis of this data. These changes 
are particularly striking in molecular and cellular biology, where incredibly complex 
interactions are revealed to be at the basis of all cell functions such as sensing, move- 
ment or reproduction. It is precisely the complexity of experimental observations 
that necessitates a more accurate and in-depth analysis. The need for a quantitative 
understanding of biological phenomena has lead to different modeling approaches. 
Many times highly advanced numerical simulations are created, and some research is 
aimed at highly realistic computer models of entire signal transduction pathways, or 
even entire cells. A very different, but equally valid approach, is to simplify possibly 
very complicated interactions to a smaller set of key components which lends itself 
easier to analytical models. 

This dissertation is concerned with models of the latter type, namely, with arriving at 
biologically meaningful results from mathematical models based on a simplification of 
experimental observations. The hope of such models is that if they do indeed capture 
the key principles of the underlying the phenomenon, then a quantitative understand- 
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ing of these principles leads to new information which has not been uncovered by the 
experiments. In successful models the mathematical analysis leads to insights which 
are unattainable (or very difficult to attain) experimentally. 

The first two chapters of this work are centered around a common theme: gradient 
sensing. Chapter El discusses a model of pattern formation due to bacteria searching 
for optimal oxygen concentrations. This work is based on experiments conducted by 
Zhulin et al. [7Sj on pattern formation of such bacteria. Our mathematical model 
presented in this chapter confirms that the experiments are produced by a novel form 
of gradient sensing, and in addition, it offers some experimentally testable predictions. 
Chapter |3] presents two models of how signal transduction events lead to the ori- 
entation of a neuron toward the appropriate target. This is an inherently difficult 
problem because of the limitations on experimental data available. In this disser- 
tation two mathematical models of neuronal gradient sensing are developed, each 
aiming at understanding a different aspect of this question. Their advantages and 
disadvantages are discussed in detail in this chapter, and it is concluded that further 
work is necessary in this area. This chapter of my dissertation is a collaboration 
with Prof. Geoffrey J. Goodhill from the Neuroscience Department of Georgetown 
University Medical Center who introduced me to the biological background of growth 
cone guidance. 

The third main topic of this work, presented in Chapter also describes sensing on 
the cellular level, but in this case the signal is not biochemical, but mechanical. The 
mathematical model describes morphological changes in the cells lining the blood ves- 
sels when they are exposed to different types of shear stresses (induced by different 
types of flow over the cell surface). The predictions of the model agree qualitatively 
with the experiments, however, improvements of the model described in this chap- 
ter are necessary in order to make quantitative predictions. This work was done in 
collaboration with John S. Tamaresis from the Graduate Group in Applied Mathe- 
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matics and Prof. Abdul Barakat from the Mechanical and Aeronautical Engineering 
Department of the University of California, Davis. 

Each chapter contains the biological terminology and data relevant to the topic, as 
well as a discussion and conclusion of the mathematical model presented. The three 
projects are quite distinct both biologically and mathematically, therefore separate 
conclusions appeared to be most appropriate. 
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Chapter 2 

Aerotaxis 

2.1 Introduction 

The study of cell motility is a broad subject with numerous applications ranging from 
understanding how nerve cells find their place in the developing brain to understand- 
ing wound healing. Of all different forms of motile cell behavior, bacterial chemotaxis 
is the best understood. Bacterial aerotaxis differs from conventional chemotaxis in 
a number of interesting ways that were highlighted by experiments conducted at the 
Loma Linda Medical School by Zhulin et al. [75] . 

Some differences between conventional chemotaxis and aerotaxis are already known, 
but the biochemical signal transduction pathways involved in aerotaxis are not. The 
purpose of this work is to demonstrate that the unusual patterns found in Zhulin's 
aerotaxis experiments are consistent with a novel form of taxis without slow adapta- 
tion. The main result is a mathematical model based on fast adaptation that charac- 
terizes aerotactic behavior and uses experimentally obtained parameters. Analytical 
and numerical results of the model are compared to experimental data. 
The biological background is introduced in Section 12.21 which explains the terms 
chemotaxis and aerotaxis in detail and describes the differences between the two. 
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Current knowledge of the biochemical pathways involved in chemotaxis is also ex- 
plained. 

A detailed description of the Zhulin aerotaxis experiments and explanation of the 
observed pattern follows, as well as questions that arise from these experiments. We 
conclude with descriptions of various chemotaxis models and their limitations when 
applied to aerotaxis experiments. In particular, we introduce the Keller-Segel model, 
and discuss the modeling assumptions. This is followed by a summary of Griinbaum's 
work on approximating a general class of equations describing random walk behaviors. 
We used his analysis to show why most conventional chemotaxis models cannot be 
applied to the Zhulin aerotaxis experiments. Some other mathematical models are 
discussed briefly (Tranquillo [75], Barkai & Leibler jlj ), and we argue that no existing 
models in the literature provide an appropriate framework for aerotaxis. 
Section 12.31 contains our mathematical model of aerotaxis. The terms of the simple 
advection-reaction equations are explained, and the main question is the determina- 
tion of the turning rates (reaction terms). A phenomenological justification is given 
for the choice of the particular terms. Appendix [X] provides the biological reasoning 
behind choosing such form for the turning frequencies and provides a simple math- 
ematical model for a receptor which could produce such turning rates. Section 12.31 
also explains how non-dimensonalization and scaling were obtained for the model. 
Section 1231 contains the results, both analytic and numerical. We show the numerical 
simulations of the aerotactic band formation. The interpretation of numerical results 
emphasizes how the model matches and predicts the band formation. Analytical 
solutions are given for the steady state of the system. Various parameter values 
which are experimentally not easily measurable are estimated. 

In Section 12.51 we summarize our findings and talks about potential future projects 
related to this topic. 
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2.2 Background 

2.2.1 Conventional chemotaxis 

Bacterial chemotaxis is a term used for motility in the direction of higher nutrient con- 
centrations (such as sugars and amino acids) and away from repellents. In a neutral 
environment (i.e. one with uniform chemical concentrations) bacteria swim smoothly 
in a given direction for a period of time, then go through a period of abrupt changes 
of direction, called tumbling. The sequence of 'runs' and 'tumbles' results in a ran- 
dom walk. This random walk becomes biased when attractant (or repellent) is added 
(or removed). Additional attractants suppress the frequency with which tumbling 
occurs; therefore, the straight runs lengthen in the direction of the highest attractant 
concentration. This allows bacteria to move up the gradient. Removing repellents 
has the same effect. On the other hand, adding repellents or removing attractants 
both increase the frequency of tumbling and facilitate the movement of bacteria down 
the gradient Figure shows how turning frequencies change in response to at- 
tractants or repellents. Keeping attractant (and repellent) concentrations constant 
results in adaptation of turning rates. 

The swimming mechanism of cells can vary greatly, resulting in different types of 
swimming not discussed in this thesis. However, the underlying principle of a biased 
random walk is prevalent in all forms of bacterial motility. 

Chemotactic movement is necessary for bacterial cells because concentration gradients 
are detected by temporal comparison. This means that rather than being able to 
measure concentrations at the different ends of the cell, bacteria must move through 
their environment to be able to detect concentration changes. In order to have a 
temporal sensing mechanism, it is crucial that cells retain some information about 
the previous environment; in other words, it is necessary for the cell to have some 
sort of memory. 
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Fast Excitation (Phosphorylation) 



Slow Adaptation (Methylation) 



add 
attractant 



add 
attractant 



remove 
attractant 



tumbling 






Figure 2.1: Changes in turning frequencies as a response to attractants and repellents. Fast excita- 
tion and slow adaptation are demonstrated. (Figure based on Bray, \l'2\) 

The mechanism for temporal sensing is also dependent on adaptation to the current 
level of attractants and nutrients. This allows cells to remain sensitive to concen- 
tration changes in a wide range of chemical environments. Adaptation and memory 
are related concepts, since memory is a consequence of a slow adaptation mechanism 
which allows cells to retain information about the previous environment for a period 
of time. During adaptation, the cell remains in a chemical state determined by the 
attractant (repellent) concentration before, and this time lag between the current 
and the past states serves as the memory of the cell. Any model of conventional 
chemotaxis must address the issues of sensitivity and memory. 

The signal transduction pathways in bacteria such as Escherichia coli have been 
widely studied jl]. There are two important chemical reactions, phosphorylation and 
methylation, which are responsible for controlling the tumbling frequency, and which 
act on very different time scales. Phosphorylation is a very fast reaction (on the order 
of milliseconds [12]) that causes a fast response (on the order of 0.1 sec) to changes 
in the attractant or repellent concentration, while methylation acts on the same time 
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scale as a single run (1-3 sec), and it is responsible for the adaptation to current 
chemical concentrations. 



attractants and repellents 




flagellum 



transport protein 



T a Pj chemotaxis receptor 



Che protein 



Figure 2.2: Signal transduction pathway in bacteria. (Figure based on Bray, |12j.) 



A simple model of the signal transduction pathway shown in Figure 12.21 involves a 
number of proteins, CheA, CheW, CheY, etc. The signaling works as follows: ligands 
(attractants and repellents) bind to transmembrane chemoreceptors of the bacteria. 
The cytosolic sides of the receptors are attached to kinases CheA and CheW. CheA 
phosphorylates both itself and CheY, which is a protein that acts as a messenger inside 
the cell. The role of CheY is to increase the tumbling frequency by reacting with the 
motor. CheY eventually dephosphorylates with the help of CheZ. Attractant binding 
slows down CheY phosphorylation, which in turn leads to suppressed tumbling and 
longer runs in the direction of attractants. Methylation, on the other hand, speeds up 
CheY phosphorylation which leads to an eventual return to the base tumbling rates, 
i.e. adaptation. Methylation occurs through a pair of enzymes, CheR and CheB, 
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which add and remove methyl groups. As the rate of methylation of the receptors 
slowly increases, the tumbling frequency increases as well, and the cell returns to 
the original turning rates. There is also a coupling between the phosphorylation and 
methylation pathways, and it is the phosphorylation of CheB by CheA. This results 
in an increase in the demethylation activity of CheB, so this also contributes to the 
suppression of tumbling rates |3|. 

2.2.2 Aerotaxis 

Aerotaxis is a specific type of taxis (directed movement) in which both the attractant 
and the repellent are particular concentrations of oxygen. Very low and very high 
concentrations both act as repellents whereas some intermediate concentrations at- 
tract bacteria. The particular range of desirable concentration depends on the species. 
Aerotaxis was first discovered in 1676 by van Leeuwenhoek who noticed aggregation 
of cells underneath the surface of a solution. Later, in 1881, Englemann observed 
that bacteria aggregate at the edges of coverglass and around the air bubbles trapped 
underneath [7T] . Initially, aerotaxis was considered a chemotactic response to oxygen, 
but further studies, summarized in Taylor's 1983 paper jTTJ , reveal several aspects in 
which aerotaxis and conventional chemotaxis differ drastically. 

First of all, we list these differences between conventional chemotaxis and aerotaxis, 
then explain their meaning and significance below. One of the most notable differ- 
ences is that aerotaxis is metabolism-dependent. This means that while bacteria do 
consume the oxygen, chemotactic bacteria can be attracted to nutrients which they 
are unable to metabolize pQ . There is also quite a bit of evidence jSUl EH EE] , that 
aerotaxis is an example of so called "energy taxis" in which cells monitor their in- 
ternal energy balance and react to optimize it. (In contrast, chemotaxis is based on 
monitoring and optimizing the nutrient availability in the external environment.) We 
discuss the notion of energy taxis further below. The third factor distinguishing aero- 



CHAPTER 2. 



MATHEMATICAL MODELS IN BIOLOGY 



10 



taxis and chemotaxis is the signal transduction pathway. This includes differences 
in the receptors utilized, as well as the fact that aerotactic movement does not have 
a methylation-dependent adaptation [71]. The precise mechanism of adaptation in 
aerotactic bacteria is currently unknown. 

In order to explain energy taxis, we must introduce some new terms. Proton motive 
force refers to the electrochemical potential difference across the membrane j7J| , and 
it is produced by linking electron transport due to respiration to translocation of 
protons [7T]. The coupling between the proton motive force and the electron transport 
system is tight, and currently it is not known which of the two acts as a signal 
for the cell's behavior ]^\. However, it is believed |oTH 173] that aerotactic bacteria 
respond to internal changes in the proton motive force or the elector transport system 
and not directly to the extracellular oxygen levels. This is supported by a series 
of experiments [7T| 173] in which other signals that changed electron transport also 
resulted in behavioral responses. For example, Taylor and Zhulin [73] cite cases in 
which metabolized substrates elicit tactic response, as do chemicals which are able to 
donate electrons to or accept electrons from the electron transport system. 




Figure 2.3: Aer and Tsr receptors. (Figure from Taylor and Zhulin, 73 ) 
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There are two receptors that have been proven to act as signal transducers . The 
two receptors, Aer and Tsr are shown in Figure 12.31 Aer is a novel receptor which 
plays no role in conventional chemotaxis. Rebbapragada et al. jSHl demonstrate that 
when mutants with a deactivated aer gene are exposed to an oxygen gradient, they 
find optimal oxygen concentrations much slower and less efficiently than wild type 
(non-mutant) bacteria. Meanwhile, the same mutants still exhibit normal chemotaxis 
indicating that the signal transduction pathways for chemotaxis and aerotaxis differ. 
When expression of Aer is restored, the aerotactic behavior returns. 
Tsr, the other aerotaxis receptor, also works as a receptor for conventional chemotaxis, 
sensing external environments as well as monitoring internal pH. The aer tsr double 
mutants were not capable of aerotactic sensing; however, upon restoring one or both 
of the Tsr and Aer receptors, aerotactic response returns The signal transduction 
pathway of aerotaxis is not well understood, although it is believed to converge with 
the phosphorylation pathway of conventional chemotaxis jHOl- CheA, CheW and 
CheY are part of the signal transduction pathway for aerotaxis (HO], but there is 
evidence that adaptation is methylation- independent [71]), and it is much faster than 
the adaptation response in conventional chemotaxis |71j . 

Aerotaxis is thought to be beneficial, because finding the appropriate concentrations 
of oxygen is essential for the metabolism of some species and can, in fact, be a more 
immediate need than finding the appropriate nutrient levels [7T] . According to Taylor 
and Zhulin |7H|, the aerotactic response might prevent bacteria from getting trapped 
in anaerobic, growth-limiting environments. In their hypothesis, bacteria living in 
conditions which support growth would mostly rely on their chemotactic behavior and 
would use aerotaxis when the maintenance of optimal internal energy levels becomes 
impossible. The mechanism for aerotaxis is expected to operate on simpler principles 
than that of chemotaxis, since it is "designed" to find a well-defined range of a single 
chemical, oxygen. On the other hand, chemotaxis allows bacteria to choose between 
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different types of nutrients of possibly different concentrations and quality, as well as 
allowing adaptation to a wide range of concentrations HE] • 

2.2.3 Aerotaxis experiments with Azospirillum brasilense 




Flagellum 



Figure 2.4: Schematic figure of Azospirillum brasilense. One flagellum is attached to the ellipsoidal 
body. 

Azospirillum brasilense is a 1-2 fim nitrogen-fixing plant-associated bacterium [TH] 
with an ellipsoidal body to which one flagellum is attached. (See Figure 12.41 for a 
schematic diagram of A. brasilense.) Counterclockwise rotation of the flagellum pro- 
duces forward motion, while clockwise rotation reverses the direction. The essentially 
one dimensional movement of A. brasilense makes it a simple organism to model. It 
is generally accepted that its positive aerotaxis (attraction) is a response to changes 
in the proton motive force, and in Zhulin's hypothesis [7H] this is also the signal for 
negative aerotaxis. 

A. brasilense is aerobic, but it prefers very low concentrations of dissolved oxygen. 
Zhulin's experiments demonstrated by spatial and temporal assays that oxygen indeed 
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acts both as repellent and attractant, and Zhulin has provided evidence j7S] that both 
are linked to monitoring the proton motive force inside the cell. As before, one can 
reason that aerotaxis provides an additional advantage for A. brasilense by guiding 
it to an optimal range of oxygen concentration. Nitrogen fixing can be accomplished 
only in environments where the oxygen concentration is below 1 percent and in 
these environments A. brasilense remains capable of maintaining aerobic metabolism 

EB|. 




30 60 90 120 25 

Time(s) Time(rain) 



Figure 2.5: The graph on the left shows the turning frequency as a function of time. Significant 
changes occur as the oxygen concentration jumps to 0.5 percent from no oxygen, and as the oxygen 
concentration jumps to 21 percent. The figure on the right shows the proton motive force as a 
function of time. (Figure based on Zhulin. [75]') 

In the temporal assays for aerotaxis, a small droplet of bacterial suspension was spread 
on a slide [7S] then exposed to different oxygen concentrations. Figure l2~5l illustrates 
the results of these experiments. It shows that the bacterial turning frequency re- 
mained constant at approximately 0.28 s^ 1 for most changes in oxygen concentration. 
There were two instances where significant changes were observed in the turning fre- 
quency. The reversal frequency dropped to 0.09 s~ l when cells were ventilated with 
0.5 percent oxygen after being equillibrated to nitrogen (no oxygen). The other no- 
ticeable change occurred when cells adapted to 0.5 percent oxygen were exposed to 
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5 percent oxygen. In this case, the reversal frequency jumped to 0.49 s 1 . Table I2~T1 
summarizes the results of these experiments. 



Change in oxygen concentration (percentage) 


Reversal frequency (s 1 ) 


21 to 100 


0.32 ± 0.04 


to 21 


0.27 ± 0.04 


100 to 21 


0. 24 ± 0.06 


21 to 


0. 22 ± 0.05 


to 0.5 


0.09 ± 0.03 


0.5 to 5 


0.49 ± 0.09 



Table 2.1: Turning frequencies as the oxygen concentration changes. (Table from Zhulin, |78j1 

These findings indicate that cells are attracted to 0.5 percent oxygen concentration, 
since both increases and decreases in the oxygen concentration caused negative aero- 
tactic response. Zhulin also calculated the proton motive force for various oxygen 
concentrations and found that the highest values of proton motive force corresponded 
to the 0.3-0.5 percent oxygen concentration range. 

In the spatial assay for aerotaxis, a flat 50 by 2 by 0.1 mm capillary tube was filled with 
a solution containing A. brasilense distributed uniformly, with no dissolved oxygen 
in the solution. At the open end of the capillary, the oxygen concentration was 
maintained at various fixed levels, and this induced a band formation inside the 
capillary within 50 seconds to 3 minutes, depending on the oxygen concentration. 
The oxygen diffused into the capillary and was consumed at a certain rate by the 
bacteria. The bacteria aggregated to the favorable oxygen concentration, and the 
region where their density was high was seen as the band inside the capillary. If 
the oxygen concentration was 100 percent, the band moved further away from the 
meniscus. If nitrogen replaced oxygen, then the band moved closer to the open end, 
and eventually disappeared. If air was introduced again, then the band reappeared. 
An unusual feature of these spatial assays was the steepness of the gradients produced. 
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Oxygen concentrations change from 20 percent to zero oxygen in about 1.6 mm. 





5 



Figure 2.6: Figure A shows the oxygen concentration as a function of space. Figure B shows the 
aerotactic band on the same scale. The band forms about 1.6 mm form the meniscus where the 
oxygen concentration is 0.3-0.5 percent. (Figure based on Zhulin, |75]1 

In the spatial assay, the swimming speed of bacteria was observed. Inside the band, 
the cells swam with an average speed of 49^ whereas outside the band the average 
speed was 14 — 22^, depending on the oxygen concentration. The cells swimming 
in either direction would swim straight through the band then reverse direction im- 
mediately once they passed the band. The bacteria did not reverse their direction 
of swimming while inside the band [7H]. The bacterial density inside the band was 
nearly a hundred times that in front of the band (between the band and the meniscus). 
The bacterial density behind the band remained approximately constant. Figure l2~6l 
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shows the bacterial band which evolves in the spatial assay, and, on the same scale, 
the oxygen concentration corresponding to the various bacterial densities. 
Our mathematical model will have to answer the following questions that arise based 
on the experiments: Can we create a model which leads to the evolution an aerotactic 
band with similar characteristics? We would like to test whether a model based on 
the simple turning rules which are suggested by the experiments would be able to 
produce the observed pattern. In particular, the hope is that of the model would 
demonstrate that the measurable quantities, such as the bacterial density inside the 
band and outside the band, the width of the band, and its distance from the meniscus, 
can be explained by assuming our simple turning rule. 

Some deeper questions to pose beyond this model would be: What kind of signal 
transduction mechanism do aerotactic bacteria have? In particular, what sort of 
receptors make it possible for oxygen to act both as an attractant and a repellent? 
What sort of properties of this signal transduction mechanism allow the cells to 
react to such steep concentration gradients? How can the turning frequency of the 
cells change so abruptly at the boundaries of the aerotactic band? What is the 
underlying mechanism that allows the turning frequency to respond immediately? 
These questions are likely to be answered by further experiments. 

2.2.4 Mathematical models for conventional chemotaxis 

This section provides a brief overview of the mathematical tools used to describe 
conventional chemotaxis and briefly discusses why these existing models cannot be 
employed to study aerotaxis. 

The first mathematical models for bacterial chemotaxis were created in the early 
1970's by Keller and Segel [1111221 • m their work [38J, they show how local gradient 
detection by individual cells produces fluctuations in their path, and how the average 
over many cells corresponds to a macroscopic flux. They derive the macroscopic flux 
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equations based on individual bacteria. Although in this paper they make incorrect 
assumptions about the way bacteria are able to detect gradients (namely, they assume 
that bacteria can compare concentrations at their head and tail), they arrive at a very 
useful expression for the macroscopic flux: 

t/ \ fdb s .dL 

J(X) = -M-7-) +X0-7- 

In this equation fJ,(L) is the "diffusion" , or motility coefficient, and x{L) is the chemo- 
tactic coefficient. b(x,t) is the bacterial density, and L(x,t) is the concentration of 
the ligands. The rate of change of the bacterial density b(x, t) can be given by 

db 

5 = "™ 

with the appropriate boundary condition. In one dimension, this becomes the well 
known Keller-Segel chemotaxis equation: 

^ - ® ( fi^ b +xb^ L ) (2 1) 

dt dx dx dx 



Again, boundary conditions must be imposed. It is clear in the equation ()2.1|) that 
the first term on the right-hand side is diffusion due to random motility, and the 
second term is due to chemotactic flux. Usually it is also assumed that the spatial 
gradient of the concentration is also small, and therefore ^ can be approximated by 
a constant, and absorbed in x, the chemotactic coefficient. 

In a later paper [S3], Segel derives the same chemotactic equation (j2.1|) based on 
changes in receptor configuration as a result of attractant (or repellent) binding. In 
this model, there is no need to use the incorrect hypothesis that bacteria are able 
to compare concentrations at different parts of their body. Instead, this model is 
developed by writing down a separate equation for left- and right-moving bacteria in 
different receptor configurations. 

In order to to arrive at ()2.1|) from this system of equations, Segel must use the 
assumption that spatial gradients are small. Almost the same analysis is summarized 
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more lucidly by Griinbaum; therefore, in order to understand where the assumption 
of small spatial gradients is needed, it is appropriate to look at Griinbaum's article 
"Advection-diffusion equations for internal state-meditated random walks" [23] . 
Griinbaum gives a simple argument regarding the validity of parabolic advection- 
diffusion approximations to hyperbolic advection equations. He starts with a system 
of one-dimensional advection equations describing left- and right-moving individuals, 
similar to the Segel system in [S3] . 

<%+ db + _ T _ +l+ 

— + v— = a b -a + b + 
at ox 

db~ db~ ,„ . 

Boundary conditions are omitted again for the time being. The turning rates are 
denoted by a + , <j~ and velocity by v. The velocity is chosen to be positive for the left- 
moving bacteria, and negative for the right-moving bacteria. Total bacterial density 
is given by b(x,t) = b + (x,t) + b~(x,t); in other words, the total bacterial density 
is the sum of the right- and left-moving terms. Griinbaum also defines the density 
flux, J(x,t) as J(x,t) = v(b + (x,t) — b~(x,t)), and introduces two new variables, 
cr = |(o~ + + o" - ) and Act = — <t _ ). By taking the sum and difference of the 

two equations in (J2.2)) and making the appropriate substitutions, he arrives at the 
following two equations: 

d ±-- d -l (23) 
dt~ dx {Z - 6) 

1 dJ T 1 . 9 db . . . 

J = - — (-v - — h 2Aavb) (2.4) 



2a dt 2a dx 

In equation (|2.4|) . the objective is to estimate the term J t using the conservation 
equation (|2.3|) . 

Griinbaum does this the following way. Assume that spatial derivatives are small, 
then a small J x implies that 6 t is small, or the bacterial density varies on a slow 
time scale. The solution to ()2.4j) is the linear combination of the homogeneous and 
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inhomogeneous solutions. The inhomogeneous terms all contain b(x, t) or a spatial 
derivative of b; therefore, they vary on a slow time scale. However, the solution to 
the homogeneous equation 

1 9J J n 

2^ + J = ° 

will be an exponential term acting on a faster time scale. This implies that after 
a short initial period, J approaches a quasi-equilibrium. This gives the following 
equation: 

J « —(~v 2 — + 2A(rvb) 2.5 
2(T ox 

By substituting this approximation for the flux into ()2.3|) . one can recover the Keller- 
Segel equation: 

db d db 

with 

- f_ 

X = (2.7) 

The advantage of this formulation over the Keller-Segel equation is clear, since ()2.7|) 
shows that the two important coefficients, the motility coefficient and the chemotactic 
coefficient, can be expressed in terms of empirically observable quantities, namely the 
cell velocity and the turning rates. It is also clear from Grunbaum's analysis that the 
assumption of a small spatial gradient is necessary in order to arrive at the Keller- 
Segel chemotaxis equation, ()2.6|) . 

In addition to the assumption of a small spatial gradient, another difficulty with 
this model is that an exact measure of the difference of turning frequencies, Act, is 
very difficult to obtain empirically, or to approximate analytically. An analytical 
expression for \i the taxis coefficient, can be derived in terms of the characteristic 
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time scale of a run, the bacterial speed and the characteristic attractant concentration. 
This analytical result does not rely on having to measure the turning frequency, but is 
rather obtained by a perturbation method which separates time scales of the reaction. 
In this perturbation analysis, the characteristic run time is obtained based on the 
assumption that adaptation is slow. This is not a correct assumption in the aerotaxis 
experiments. 

Another significant development in mathematical modeling of chemotaxis came in the 
late 1980's and early 1990's. A large effort was made to unite knowledge about recep- 
tor dynamics and signal transduction pathways (many times referred to as 'internal 
state dynamics') to random motility fHIO]- I n many of these models (e.g. |51|). 
not only bacterial chemotaxis but chemotaxis of animal cells is considered, making 
the internal state dynamics far more complex than in the bacterial case. There were 
several models, for example those by Barkai and Leibler 0] and by Tranquillo et al. 
|18| I5T] which included very detailed models for the biochemical mechanism. 
One such detailed model of chemotaxis in an animal cell is due to Tranquillo [SI]. The 
general approach of Tranquillo and his collaborators is as follows. Identify a simpli- 
fied (but still realistic) scheme to describe all possible receptor states and intracellular 
chemicals that govern the turning behavior. All rate constants and parameters are 
approximated based on empirical data. Each of the receptor states are interdepen- 
dent stochastic variables, and their time evolution is calculated using a multivariate 
probability density function. 

The jump processes between the various states can be approximated by continuous 
stochastic differential equations by making certain assumptions about the system. 
The probability function must be linearized in order to be represented by an analo- 
gous system of stochastic differential equations. However, the linearization procedure 
implicitly assumes small fluctuations in the stochastic variables. This, again, is an 
assumption that would be violated in the aerotaxis experiments. 
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In order to analyze the stochastic differential equations, deterministic equations are 
derived either using averages of the stochastic variables or transformations from the 
stochastic variables to new, deterministic variables. With this method, the Fokker- 
Planck equations are derived and used to estimate cell movement on longer time scales 

EH- 

These models ([ISIEI]) are highly sophisticated biologically and mathematically alike, 
although lack the satisfying simplicity of the earlier Keller-Segel equations. As men- 
tioned above, they also assume certain features of the internal state dynamics that 
would be violated in steep attractant gradient; therefore, they cannot be used for our 
purposes. 

Barkai & Leibler also create a model which is very closely built on experimentally 
obtained data. In "Robustness in simple biochemical networks" @] they examine 
the question of receptor adaptation in bacteria, and they propose a quantitative 
model in which a wide range of biochemical parameters are admissible. The model 
Barkai & Leibler propose is a system of ordinary differential equations that is based 
on the accurate description of the possible receptor states. The key to the model 
is a small network that has two states: an active and an inactive state. In the 
active state the external signal leads to a fast response, whereas in the inactive state 
there is no response. The shift between the active and inactive states is a slower 
process, which corresponds to the receptor methylation. The system exhibits perfect 
adaptation if the active state is independent of the magnitude of the outside stimulus. 
A model based on the same principles is presented in the chapter on animal cell 
chemotaxis, Section 13.3.21 The contribution of Barkai & Leibler is significant for 
bacterial chemotaxis, because their work shows that fine-tuning the model parameters 
is unnecessary in order to achieve adaptation. However, in our model for aerotaxis 
there is no methylation adaptation, and we argue that if there is adaptation, it must 
occur on the fast time scale. 
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There is also some literature that focuses directly on aerotactic behavior and models 
for aerotaxis [201 EH]. I n [HI]; certain bioconvection patterns are described quanti- 
tatively by developing a model of aerotaxis. Kessler's experiments described in this 
paper involve an initially well-stirred suspension in which cells swim upwards toward 
oxygen, then, after the top layer becomes sufficiently denser than the bottom layer, 
an instability occurs. The overturning instability evolves into the observed patterns. 
Some of the phenomenon is similar to the Zhulin aerotaxis experiments, namely, here 
bacteria consuming oxygen create the gradient while the oxygen level is fixed at the 
meniscus. However, in the Kessler experiments the bacterial convection stirs the 
solution, and the bacteria carry oxygen into deeper layers of the solution. 
The mathematical model describing these experiments consists of a conservation equa- 
tion for the cells and a reaction-diffusion equation for the oxygen concentration. The 
Navier-Stokes equations are not required, since no bulk fluid flow is assumed. The 
authors arrive at the same equation as the Keller-Segel equation for the bacteria. The 
whole system is: 



In this model L(x, y, z, t) is the concentration of oxygen, b(x, y, z, t) is the density of 
bacteria, as above, u is the fluid velocity taken to be zero, v(Q) is the average cell 
velocity which is a function of a dimensionless measure of L denoted 0, D and Dl are 
the diffusion coefficients of the bacteria and of oxygen, respectively, and, finally, k(Q) 
is the coefficient of nutrient consumption by the bacteria. The boundary conditions 
are applied at z = and at z — — h, the bottom of the suspension. They are given 





dL 
~dt 



V ■ (Lu - D L VL) - kb 



(2.9) 
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by: 



L = L at z = 0, — — = at z = —h 

oz 

db 

v z b — D— = at z = and —h 



(2.10) 



(2.11) 



These boundary conditions mean that there is no flux of oxygen or bacteria at any of 
the boundaries. The initial condition is given for a well-stirred solution and uniform 
suspension: 



By non-dimensionalizing the equations, the authors arrive at a perturbation problem 
whose analysis leads to a good quantitative description of the experiments. There are 
several reasons why this model would not be valid for the Zhulin aerotaxis experi- 
ments. First, just like the Keller-Segel equation, Hillesdon et al. make use of the small 
gradient assumption implicitly. Also, the Kessler experiments are examples of kinesis 
rather than taxis. In kinesis the movement is determined by local concentration of 
attractant, not by concentration gradients, as it is in taxis. 

In his review article [30] . Hill summarizes other models of various types of chemotaxis 
related to pattern formation (gyrotaxis, geotaxis, phototaxis). These equations make 
very similar assumptions to the models discussed above, namely, conservation equa- 
tions for the bacteria involving gradients of the flux due to swimming and random 
motility. When, in addition, randomness in drift speed and direction of motility are 
introduced, the Fokker-Planck equation is used to describe tactic behavior. 
In all the above mentioned models, there are two important aspect of chemotaxis that 
must be taken into account: the gradient of attractant concentration and the adapta- 
tion time. For small gradients and fast adaptation, it is simple to find approximations 
to the hyperbolic advection equations describing bacterial motion. Grunbaum also 




(2.12) 



b(z, 0) = b 



(2.13) 
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shows that for slow adaptation in shallow gradients, one is able to simplify the equa- 
tions to a parabolic system. However, if the attractant gradients are large, then these 
approximations do not work for slow adaptation any more. In this case, the only pos- 
sible approach is to try Monte Carlo simulations. The model of aerotaxis presented 
in the next chapter uses the fact that while gradients are large, the adaptation time 
in our case is fast. This allows us to solve the original hyperbolic system. 
It is also possible to give a heuristic argument for why conventional chemotaxis models 
that rely on slow adaptation cannot be used to model the Zhulin aerotaxis experi- 
ments. The following Monte-Carlo simulation-Figure 12.71 illustrates this reasoning. 

Dashed line: path of a bacteria, solid line: turning frequency 

7 

6 

5 

4 

3 

2< 

1 



-1 

-2< 
( 

Figure 2.7: Monte-Carlo simulation. 

The wide solid lines (at -2 and 2) represent the sides of the capillary tube (so the 
bacteria are confined to this region), and the dotted lines (at -1 and 1) represent the 
favorable oxygen concentrations. The following rules govern the movement of each 
cell: 

• a cell moves straight to the left or to the right with a constant velocity, v; 

• turning frequency inside the favorable region (between -1 and 1) is a = 0; 
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• the cell leaves the band at some random time, r, and at this time the turning 
frequency jumps form a = to er = c; 

• characteristic time of adaptation is t a ; 

• adaptation to the baseline turning frequency of a = is exponential, and it is 
given by a = ce *» . 

The adaptation in these simulations is assumed to be slow, which means that |, the 
time of a straight run, is of the same order of magnitude as t a , the characteristic time 
of adaptation. Later, in our model of aerotaxis we assume that t a « K 
In these simulations, the particle's velocity and turning frequency are given deter- 
ministically. When the cell is outside the band, the time of turning, r, is determined 
based on the difference of a uniformly generated random number between and 1 
and the turning frequency. 

In the figure, we can see a typical run of the simulation. Most of the time the bac- 
terium stays inside the optimal oxygen concentration, because upon leaving the band, 
its turning frequency jumps from to c, and it is likely to turn back into the band. 
However, outside the band the turning frequency is large for a period of time (t a ) 
due to slow adaptation, and it frequently causes the cell to keep tumbling and get- 
ting trapped outside the optimal environment. Running the simulation up to 10,000 
times, the bacterial density inside the band was only three times the density outside 
the band. This is clearly very different from the 100:1 ratio observed experimentally. 
One must conclude that there are no existing models of bacterial chemotaxis (other 
than Monte Carlo simulations) that can describe the behavior in steep gradients. Most 
existing chemotaxis models also assume a slow adaptation of the turning rates. Exact 
mathematical descriptions of turning rates based on slow adaptation are very difficult 
to analyze, and in order to create tractable equations, one must make approximations. 
The approximations involve assumptions of small spatial gradients, since this allows 
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continuity of internal state variables. Assuming slow adaptation and a steep spatial 
gradient, no approximations are possible leading to simple mathematical expression. 
This suggests that it would be futile to attempt to model the Zhulin experiments with 
already existing chemotaxis equations. However, since aerotaxis is known to have fast 
adaptation, mathematical expression of the turning rates is much simpler; thus, we 
can develop a different model in which one need not rely on approximations based on 
shallow gradients. 
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2.3 Model 



2.3.1 Mathematical model for aerotaxis 

Now we can present the mathematical model for aerotaxis which describes the pattern 
formation as a result of steep concentration gradients. A brief derivation of the 
equations is given below. The bacteria's movement is governed by advection-reaction 
equations. The advection term describes the directed swimming of bacteria, while 
the reaction terms denote the turning of bacteria in response to the oxygen gradients. 
The evolution of the ligand (oxygen) concentration is given by a reaction-diffusion 
equation which is coupled to the advection equations through the term describing 
oxygen consumption by the bacteria. 

We assume that bacterial movement is one-dimensional. Although some turning is 
possible, this assumption is based on empirical evidence of A. brasilense swimming. 
The cells are observed to swim either forward, or, upon changing their direction, back- 
ward. The hypothesis that the turning rates are dependent on oxygen concentration 
and the oxygen gradient is the main assumption of the model, and it is discussed at 
length below. The boundary conditions imposed simply mean that all left-moving 
cells turn to the right at the left boundary, and similarly, all right-moving cells turn 
to the left at the right boundary. The assumption of the conservation of the num- 
ber of bacteria is used later, and it justifies the lack of birth and death terms in 
the equations. Since the band forms on the order of minutes, this is a reasonable 
omission. 



dr 
di 




dx 



~ frlT + flrl 



(2.14) 



dl _ d(vl) 



+ frir - flrl 



(2.15) 



dt dx 




d 2 L 



k(r + /) 



(2.16) 



dx 2 



CHAPTER 2. 



MATHEMATICAL MODELS IN BIOLOGY 



28 



r(0) = 1(0) 
r(a) = 1(a) 

r(x,t) - right moving bacteria 

l(x,t) - left moving bacteria 

L(x,t) - ligand (oxygen) concentration 

fri(L) - rate of turning from right-moving to left-moving cell 

fir(L) - rate of turning from left-moving to right-moving cell 

v(L(x)) - bacterial speed 

D - diffusion coefficient for oxygen 

k - rate of oxygen consumption by bacteria 

a - length of the capillary 

The initial condition for the left- and right-moving bacterial densities is the same 
constant for all positions, x, such that the sum of the two populations is the total 
bacterial density. The initial condition for the oxygen concentration is L Q at the left 
boundary, and zero everywhere else. 




x x+A x 

Figure 2.8: A volume element that bacteria swim through. 



A brief justification (similar to the reasoning of Segel in 63J) of using the advection 
equation (|2.14j) is as follows. Consider the cells swimming in a capillary of a fixed 
cross-sectional area, R. Let us look at a short section of the capillary, from x to 
x + Ax. The density of right-moving cells in this section is given by rRAx. Then 
the rate of change of the r, the right-moving bacterial density, with time respect to 



CHAPTER 2. 



MATHEMATICAL MODELS IN BIOLOGY 



29 



time is equal to (i) change due to reversal of direction to become left-moving, and, 
similarly, left-moving bacteria turning to become right-moving; (ii) change due to 
cells swimming into and out of the slice. 

The change due to turning is quite straight-forward. It is the turning rate of cells 
times the bacteria in the given volume, or fi r lRAx — f r irRAx. The term due to cell 
swimming is the rate at which cells flow into the cell, and the rate at which they flow 
out, 

v(x)r(x)R — v{x + Ax)r(x + Ax)R. 

This leads to 

dirRAx) 

= fi r lRAx — f r irRAx — R[r(x + Ax)v(x + Ax) — r(x)v(x)}. 

Dividing by RAx and letting Ax approach zero, we obtain 

dr d(vr) 

^7 = « Jnr + Iirl 

at ox 

Equation (|2.15|) can be obtained in a similar fashion. 

The most important question regarding the model is the determination of the the 
turning frequencies. There are two questions to be answered: what function of L are 
the turning frequencies? What is the biological evidence in support of this mathe- 
matical expression? Both of these questions are answered below. 

2.3.2 Expression for the turning frequencies 

We must consider the mathematical expression for the turning rates. It is known 
from the experiments that there is a range of concentrations where bacteria do not 
turn around. We can call this range L min to L max . This is the range where the proton 
motive force (PMF) is the highest; therefore, it is the preferred concentration range 
where the bacterial band develops. Outside this range, there must be some threshold 
values, L min and L max such that if a cell is between L min and L min , or between L max 
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and L max , it turns back inside the band. The bacterial frequencies are shown in 
Figure El 




Figure 2.9: Turning rates of the bacteria. 

However, the turning rates cannot depend on the values of oxygen concentration 
alone, since if that were the cell getting to L min from outside the band would 

have to turn, and would therefore never enter the favorable oxygen range. This 
suggests that the bacteria must be able to retain some additional information about 
the environment, for example the gradient of the ligand. This would allow a cell 
arriving to L m i n to determine whether to keep on swimming (if it came from outside 
of the band) or to turn around (if it came from within the band). 
In our hypothesis, this additional information comes from monitoring the proton 
motive force inside the cell. As it is shown in Figure 12.101 the PMF has its largest 
value, and it is a constant, when the cell is between L min and L max . Outside the 
band, the PMF has a low value, and it is also a constant. The turning signal for a cell 
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min 25 

Time(min) 



Figure 2.10: The figure shows the proton motive force (PMF) versus the oxygen concentration. The 
highest proton motive force is observed between L m i n and L max . PMF is increasing between L m in 
and L m i n and decreases between L max and L max . (Figure based on Zhulin, |78j.) 

is the negative temporal gradient of the PMF. Positive temporal gradients and the 
high constant value of PMF suppresses tumbling. Cells are able to detect temporal 
gradients of PMF by swimming through a spatial gradient of oxygen (which is linearly 
proportional to time). The biological justification for such turning rates is given in 
Appendix [X] Further details on a model of a receptor producing such turning rates 
are also included. 

The easiest way to construct such turning rates is by defining piecewise linear func- 
tions as is shown in the above figure (Fig. 12. 9|) . We express the turning rates of the 
left-turning bacteria and the turning rate of the right-turning bacteria separately, so 
now the two rates only depend explicitly on the ligand concentration, and their de- 
pendence on the oxygen gradient is implicit. (If a general turning rate was given for 
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all bacteria, this rate would explicitly depend on the gradient of oxygen.) This simple 
choice for the turning frequencies makes the system of equations almost linear, with 
the only non-linearity resulting from the dependence of turning rates on the ligand 
gradient. Now we can explicitly write down the turning rates. 



rl 



fl 



Ir 
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L *C L m i n 
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c, 
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c, 


Lmax *^ 



(2.17) 



(2.18) 



In all the above equations C is some larger turning rate than c. The exact value of 
C and c will be specified later. Now with equations (ETHJ) . (ETToT) . (EHTfl) . (l2~T7j) and 
([2.18)1 . we have the complete system of equations describing bacterial swimming. 



2.3.3 Non-dimensionalization and scaling 

We have the full system of equations for the aerotaxis experiments including turning 
rates and boundary conditions. 

dr d(—vr) 
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c, 
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r(0) 


= /(0) 




r(a) 


= J(a) 



Before the system can be analyzed, it is important to non-dimensionalize all variables. 
The following parameters are important: 



Size of capillary tube: 50x2x0.1 mm 
Preferred [0 2 ] : 0.3-0.5 percent 
Band width: 0.2 mm 

Distance of band from capillary tube end: 1.6 mm 
Time of band formation: 50 sec - 3 min 
[O2] outside the capillary tube: 21 percent 
Speed: 40^ 

1 sec 

Diffusion coefficient: 2 • 10" 9 ^ 

sec 

Turning frequency: lsec -1 

Rate of oxygen consumption: 3 • 10~ n ( ee ^ ee ) 
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The appropriate spatial scale can be found by estimating the distance which would 
allow the bacteria to outrun the invasion of the oxygen. The diffusion of oxygen 
implies that x ~ Voi, whereas the distance for the escape of bacteria is x ~ vt. From 
these two expressions we can approximate the time for the cells to outrun diffusion, it 
is to ~ ~ 5 sec, which means that Xq ~ v to ~ 100 /im. This suggests that the time 
scale should be on the order of about 10 seconds (which agrees with the experiments 
where the band develops between 50 seconds and 3 minutes). For the spatial scale, 
we choose 2 mm, since this is the length of the region where the bacteria are found, 
but from the scaling argument it is clear that the resolution must be smaller than 100 
fim. Based on this we arrive at the characteristic scales summarized in Table 12.21 



Measurement 


One unit 


Length 


2 mm 


Time 


10 sec 


Oxygen concentration 


i m. 

mi 


Bacterial concentration 


2-10 7 ^i 

ml 



Table 2.2: Table gives the units of measurement for length, time, oxygen concentration and bacterial 
concentration. 



Measurement 


Dimensional quantity 


Non-dimensional value 


Speed 
Diffusion coefficient 
Turning frequency 
Rate of oxygen consumption 


mum 

sec 

2 . 10 -9mi 

sec 

lsec -1 
q m-ii ^ M 

iU (cell)(sec) 


0.2 
0.01 
10 

4- 10~ 3 



Table 2.3: Non-dimensional parameter values. 



These scales give us the non-dimensionalized values for our parameters, shown in 
Table O 

Values for L min and L max are not readily measured experimentally, so several different 
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parameter values were used in the simulations. The non-dimensionalization gives the 
following differential equations: 

dr d(—r) 



dt dx 

dl _ d(l) 
dt dx 
dL d 2 L 



dt dx 2 



~ frir + flrl 
+ frlT - fi r l 

«(r + 1) 



f 



rl 



fl 



Ir 



C, L < L m in 

c 5 L m i n < L < L m i n 

c i L rn i n < L < L max 

C" T <r J <r f 

^ ; ^max \ -i^ \ ^max 

C i L max < L 

C, L < L min 

C j Lmin <^ L <^ L m in 

c i L m i n < Z/ < L max 



r> J 

j ^max 



<^ L <C L ma x 
C, Lmax < L 

r(0) = Z(0) 
r(l) = Z(l) 

Here, k = where is the rate of oxygen consumption by bacteria, t is the time 
scale, 6o an d L are the scales of the original bacterial and oxygen concentrations, 
respectively. The non-dimensional values for the turning rates, C and d are given by 
C = Ct , and d = ct . 
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2.4 Results 

2.4.1 Numerical simulations 

This chapter describes the main results of the model. The most important result is 
the series of computer simulations showing the band development as it is observed in 
the actual experiments. 

Analytical solutions for the steady state are possible, but numerical results are given 
to see the time evolution of the solutions. In order to have numerical stability, the 
equations describing the left-moving bacteria were solved with a forward-differencing 
scheme, the equations for the right-moving bacteria with a backward-differencing 
scheme. The diffusion equation for the oxygen was discretized with a forward-time, 
centered space (FTCS) scheme. The boundary condition at the right is that all right- 
moving bacteria become left-moving, and at the left boundary all the left-moving 
bacteria become right-moving. The domain was discretized by 40 grid points, repre- 
senting the full length of the capillary, 2 mm. The number of grid points was chosen 
to be 40 so the developing aerotactic band would have sufficient resolution. The size 
of the time step, At = 0.01 was chosen such that the solution to the diffusion equa- 
tion would be stable. The initial conditions in all simulations were that the bacterial 
density (both left- and right-moving) is a constant, scaled to 1. The initial oxygen 
concentration is also scaled. 

2.4.2 Numerical results 

In the following figures, one can follow the development of the aerotactic band. 
In the figures one can follow the development of the aerotactic band. Initially, (Figure 
12.11)) the bacterial density is uniform everywhere. The oxygen concentration is a 
constant at the first grid point, and zero on all other gridpoints. The apparent 
gradient of oxygen is due to the software package, Matlab connecting the first and 
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Bacteria: solid line, oxygen: dotted line, t=0 s 

2.5 I 1 1 1 1 1 



1.5 - 



0.5 - 

l 

l 

J — i , , , , , , , 1 

5 10 15 20 25 30 35 40 

Figure 2.11: Initial condition. Before the aerotactic band forms, the bacterial density is uniform. 
(Bacterial density is given by solid line.) Oxygen is zero everywhere, except at the left boundary. 
(Oxygen concentration is given by dashed line.) 



Bacteria: solid line, oxygen: dotted line, t=5 s 
3.5 | 1 1 1 1 1 
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Figure 2.12: Bacteria beginning to aggregate at the favorable oxygen concentration. Sharp oxygen 
gradient is developing. 
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Bacteria: solid line, oxygen: dotted line, t=15 s 




5 10 15 20 25 30 35 40 



Figure 2.13: After 15 seconds the band is clearly visible. Bacterial density in front of the band still 
changing. 



Bacteria: solid line, oxygen: dotted line, t=50 s 

14 1 1 1 1 1 1 
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Figure 2.14: Aerotactic band after 50 seconds. 
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Bacteria: solid line, oxygen: dotted line, t=1 80 s 
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5 10 15 20 25 30 35 40 

Figure 2.15: Quasi steady state of bacterial band. Ratio of bacterial densities in front of the band 
and inside the band agree with experimental measurements. 

second gridpoints. After 5 seconds, in Figure 12.121 the cells at the open end of the 
capillary that are exposed to high levels of oxygen concentration start swimming 
toward the lower, optimal oxygen concentrations. Meanwhile, some of the bacteria 
at the back that are close enough to the oxygen are also able to detect the optimal 
oxygen range, and the aggregation begins from both sides. Oxygen also begins to 
diffuse through the solution, and a sharp gradient evolves. 

In the figure which shows the aerotactic band after 15 seconds fFigure f2.13J) . it is 
clear that most cells from the open end of the capillary have already aggregated to 
the band. There is low bacterial density on either side of the band, because all cells at 
these positions are able to swim into the optimal range. Cells at the back of the band 
never sense the oxygen, and the bacterial density here remains practically unchanged. 
Once the band reaches its steady state (after 3 minutes), it is clear that almost all 
bacteria from the front have aggregated to the band, fFigure r2.15|) . The ratio of the 
cell density inside the band to the front of the band is more than a 100, and the 
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density inside the band to behind the band is about 10. Both of these values agree 
with the experimental data. The time scale of the band formation, the width of the 
band and the distance of the band from the meniscus also give good agreement with 
the experimental measurements. The total bacterial density is conserved in all the 

simulations. L min = 0.2 and L max = 0.7 were used in the numerical computations. 

10 i 



L m ax =4% 



6 - L max =7% 



2 -- 



L mm =0% 



0.5 



1.5 




Figure 2.16: The top figure shows how the final bacterial band changes as L max is changed. The 
bottom figure shows the numerical simulations with L max fixed, and L m m changing. 



The simulations can also help to determine the optimal values of L max and L min . 
(Figure 12.16(1 This is useful, because these parameters are difficult to measure ex- 
perimentally. By increasing the value of L max , the bacterial density inside the band 
increases too, and the the aggregation from the front of the band is more complete. 
This happens because a larger L max effectively increases the region in which cells are 
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in favorable conditions. The cells do not have to swim very far from the meniscus to 
be able to sense the favorable oxygen concentration, thus more of them get trapped 
in the optimal zone. 

Lowering the value of L m j n has a similar effect. The bacterial density is higher for low 
values of L min , and the aggregation from the back of the band is more complete in 
this case. As before, this happens because cells now get inside the band by detecting 
lower oxygen concentrations. For a higher value of L min , sensing the same oxygen 
concentration would act as a repellent, where in this case (for the low value) it becomes 
an attractant. This change in L m i n also results in the main bacterial density shifting 
toward lower oxygen concentration. 

2.4.3 Analytical results 

Here we obtain analytically explicit expressions for asymptotically stable stationary 
distributions of bacteria and oxygen. These distributions are possible to find because 
the nonlinear equations f)2. 141 12.151 12.161 12.171 12.18J) of the energy taxis model are 
piece-wise linear. In the steady state, the system of differential equations can be 
solved on intervals and replaced with algebraic equations. We solve the equations in 
two different cases. 

In the first case, we look at the fully developed band, after all the bacteria have 
aggregated here. (An estimate below shows approximately how much time must 
pass for this case to be valid.) In this case we can look at three regions, and solve 
the system of ordinary differential equations in these regions, using continuity and 
conservation arguments to match the solutions at the boundaries. In the general case, 
we assume that L max < L for the outside oxygen concentration, L , and L min = 0. 
In the next section, ( "Special steady state solutions" ) we discuss situations in which 
the above conditions do not hold. We provide solutions for L min < L < L max and 

J-'min 
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In the second case, we look at what happens if we only wait 3 minutes (in other 
words, wait just enough time for the bacteria to develop the band.) In this case, it 
is possible to show that the bacteria behind the band have not had sufficient time 
to feel the effect of the oxygen gradient; therefore, we can assume that the bacterial 
concentration there does not change at all. This allows us to solve the equations in 
two regions only, because the region behind the band does not need to be considered. 
The second case is discussed in the section "Quasi steady state solutions" . 

Steady state solutions 



0.2 0.4 0.6 0.E 



1 1.2 1.4 1.6 1.t 



Figure 2.17: Steady state distribution of the oxygen concentration and the bacterial density for the 
case when L > L max > L min . 



Let us now first consider the steady state solution when L max < L and L min = 0. 
There are three regions in this case: region I, in front of the band, from Xq = to 
x\ = d; region II, inside the band, from x\ = d to X2 = d + h\ and region III, behind 
the band, from x 2 = d + h to x% = d + h + z. Figure 0.171 illustrates this scenario. 
Let b(x) and L(x) be stationary linear densities of cells and oxygen, respectively, 
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where the x-axis has its origin at the meniscus and is directed inside the capillary. 
Let the oxygen concentration reach values L min and L max at distances X\ and X2 from 
the meniscus, respectively. Let s = v/(f r i — fi r ) ~ 100/xm, k = k/D. We use s as our 
length scale. 

In region I the bacterial density does not change with respect to time and is therefore 
given by: 

-vr' - f H r + fi r l = 0, vl' + f rt r - fi r l = 0. 

By noting that the steady state densities of left- and right-moving bacteria are the 
same and adding the two equations, we obtain: 

r = I ~ e x/s . 

Let B be the density of bacteria inside the band and d be the distance of the band 
from the open end of the capillary. Then, 

b = Be {x ~ d)/s 

The equation for oxygen in this region is: 

L" = kb = kBe {x ~ d)ls 

Here we can use the fact that the oxygen concentration outside the capillary is L Q , 
we integrate to find the oxygen concentration in this region, Lj: 

Li(x) = L - Cl x + fcBs 2 [e< x-d)/s ~ e" d/s ]. 

In region II we can write the same expression for the bacterial density as the expression 
in region I: 

-vr' - frir + f^l = 0, vl' + frir - firl = 0. 
We note that the bacterial density in this region is a constant; therefore, we have: 

r = l = B/2; b = B. 
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The equation for oxygen in this region is given by: 

L" = kB 

We integrate this expression and use the information that at one edge of the band 
the oxygen concentration is L max and the fact that the band is d distance away from 
the meniscus to find: 

L H (x) = L max - c 2 {x -d) + ^kB(x - d) 2 . 

In region III: 

vt' f r \T 4~ f\ r l = 0, vl' + f rl r-f lr l = 0; r = l~e~ x / s . 

Then, we use the continuity of the bacterial density at the boundary of the band, x 2 
to write down the equation for the density in this region: 

b = Be {x2 ~ x)/s . 

For oxygen we substitute this expression into the bacterial density: 

L" = kb = kBe {x2 ~ x)ls 

To find the equation for oxygen we integrate, and use the fact that the oxygen con- 
centration at X2 = d + h is L min : 

L IH {x) = L mm - c 3 (x -d-h) + kB S 2 [e^ d+h -^ s - 1]. 

Here we have 7 unknown parameters: B,ci,c 2 ,c 3 ,d,h,z. To find them we use the 
following boundary conditions: 



1) Continuity of bacterial density at x = x 3 : Be i{d+h) - {d+h+z))/s = e^ z ^ = b . 

2) No oxygen behind x = x 3 (which provides constant bacterial density there): 
L(d + h + z) = 0. 
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3) No flux of oxygen to x > x 3 : L'(x 3 ) = 0. 

4-5) Continuity of oxygen at the boundaries between the zones: Li{x\) = L max and 
Lii(x2) = L min . 

6-7) Continuity of flux of oxygen at the boundaries between the zones: L'j(xi) = 
L' H (xi) and L' n (x 2 ) = L' ni (x 2 ). 

These conditions lead to the following system of algebraic equations using the notation 
A = e z l s : (Note: Throughout the following calculations, we will use k = 0.003, b = 
2, L min = 0.003, L max = 0.005, s = 1.) 



B = b X (2.19) 

L m in ~ c 3 z + kBs 2 [A" 1 - 1] = (2.20) 

-c 3 + kBsX- 1 = (2.21) 

L - Cl d + kBs 2 [l - e~ d/s ] = L max (2.22) 

L m ax ~ c 2 h + —kBh 2 = L m i n (2.23) 

-a + kBs = -c 2 (2.24) 

-c 2 + kBh = -c 3 + kBs (2.25) 



From (jUlIl) and using (l2~TT^ : c 3 = kb s. From (l2~^4ll2~23|) : 

a = kb (s + hX), c 2 = kb (s + hX — sX). 

We can find analytic approximations for z, d and h. We start by approximating z. 
Substituting a, c 2, c 3 as functions of z, h into (|2.20|) . we find: 

z/s = {L mm + kb s 2 (l - e z ' s ))/(kb s 2 ) 

Let a = + 1, and ( = z/s. Then we can rewrite the expression for z/s as follows: 
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Approximating this expression we find that 



z ~ as 

(The approximation can be obtained as follows: We assume that << 1. Then we 
can express ( in the form ( = a — e, where e is a small parameter. This parameter 
can be found from: e = e~ a+e . In the zeroth approximation e ~ e~ a . Thus, z ~ as.) 
Numerically, z ~ 1.5 and e -1 » ~ 0.2231 which verifies the original assumption of 
e< « 1. 

Next, to find h we use (J2.23|) to obtain the expression: 

(v.)'-2^(v.)- 2(j x, L r ] ^ 

From above, A = e z l s ~ 4.4817, so we have ~ 1. Let y = (h/s), u = 
2{Lm kb Q ~ s *r m) - 0- 1488 - Then , from the equation 

y 2 - 2y - u = 

using the smallness of u in comparison with 1 we arrive at two roots, yi ~ |, ?/2 — 2— |. 
Only the second root corresponds to a biologically meaningful situation, since the 
width of the bacterial band must be of the same order as the length of the straight 
runs. Choosing the first solution would predict the band width to be much narrower; 
therefore, the bacteria would always have to swim across the band without turning. 
If we substitute y2 into the expression for h, then we get, in terms of the original 
variables, h = s(2 - (Lm ^~f A "" n) ) ~ 1.8512 
Finally, we approximate d/s. From (J2.22)) : 



kb Xs kb Xs z 



Since the value of d is expected to be large in comparison with s, we can approximate 
the above expression by letting e~ d l s ~ 0. We obtain: 

Lo—L m a X I \ La I \ Lp I 1 



1 + Xh/s Xh/s h/s 
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We have + 1 ~ 8.4377 >> 1 Then, using h~2s,d = 4.2188. This also verifies 

the assumption that e~ d ^ s ~ 0. 

In summary, we get the following values: 

• h ~ 2s ~ 185/im 

• z ~ (1 + j$$)s — 150/im 

m d ~ 2A f° . s ~ 370— 1860yum using the values L = 0.2 and L = 1.00, respectively 




Figure 2.18: Spatial assays for different oxygen concentrations. A, 100 percent. B, 21 percent. The 
width of the band is independent of Lq. The distance between the meniscus and the band increases 
with increasing Lq. (Figure from Zhulin, \7S\) 

Conclusion: 

(1) Because of the absence of detailed measurements of bacterial density profile, we 
do not have data on z. 

(2) The band's width, h, is close to 200/im, exactly as predicted by our theory if the 
favorable value for velocity is chosen. 

(3) The band was observed to form 400 to 900 \xm from the meniscus corresponding 
to 20 percent to 100 percent oxygen concentrations, which is in agreement with our 
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rough estimate of 370 to 1860 fxm 

(4) According to the estimate, d is proportional to L . In the experiment, d was 
observed to grow at greater values of L . Thus, we capture the qualitative dependence 
of d on L . 

(5) h does not depend on Lq. 

Both conclusions (4) and (5) are in agreement with Figure 12. 181 from Zhulins's paper, 

EB|. 



Special steady state solutions 

As discussed above, in this section we look at special cases of the steady state solution 
where the assumption that Lq > L max > L m i n no longer holds. This section extends 
the analytical results to two other cases, namely to the case were L min < L < L max 
and when L < L min . 




0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 



Figure 2.19: Special steady state distribution of the oxygen concentration and the bacterial density 
for the case when Lq < L max . 



First case: L falls between L min and L max , as shown in Figure 0.191 In this case we 
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expect to find two regions only. In region I, the band is from X\ = to x 2 = h, and 
region II, from x 2 = h to x 3 = h + z, is behind the band. 

In region I the equation for the bacterial density is a constant; therefore, it is given 
by: 

-vr' - f rl r + f lr l = 0, vl' + frir - f lr l = 0; r = l = B/2 

Here B is density of bacteria inside the band. The equation for oxygen in this region, 
using L(0) = L , and substituting the above expression for the bacterial density: 

L" = kb = kB/2, Li{x) =L + c x x + l/2kBx 2 

The other boundary condition, L(h) = L min helps to determine C\. 

L I = L + ( Lmm ~ Lo - \kBh)x + \kBx 2 

In region II we have the following equations. For the bacterial density we have 
exponential decay: 

-vr' - f rl r + f lr l = 0, vl' + f rl r - f lr l = 0; r = I ~ e"^'. 

Then, using the continuity of the bacterial density at x 2 = /i we get: 

6 = Be {h ~ x)/s 

For oxygen we have the following equations: 

L" = kb = kBe {h ~ x)/s . 

By integrating the equation twice, we have: 

L H {x) = kBs 2 e (h - x)/s + c 2 {x -h) + c 3 

Using L(h) = L min , L(h + z) = we get 

L n {x) = kBs 2 [e^ s + K 8 [ ™ (x -h) + L mrn - fcfis 2 
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To find B,h, z we need the following conditions: 

(1) Continuity of bacterial density at x — x 3 . 

(2-3) Continuity of flux at the boundaries of the regions. 

L' I (x 2 ) = L' n {x 2 ),L' n (x 3 ) = 



B = b e z/s (2.26) 
L mm -L + l kBh = _ kBs + kBs 2 (l-e-^-L mm ^ 



h 2 



r/ , kBs 2 (l - e- z / s ^ - L min 
-ksBe~ z/s + i — = (2.28) 



Substituting the expression for B into (|2.28j) we get: 

kb s 2 e z/s - kb s 2 - L min 
—kb s H = 

z 

Simplifying this equation we arrive at: 



s kb s 2 

Recalling that a = 1 + and £ = z/s we can simplify the above expression to: 

e c - C - a = 

An analytical approximation of £ is difficult to achieve, since all the terms of the 
equation are expected to be of the same order of magnitude. Based on a numerical 
estimate, ( = 0.85. This, indeed, verifies that all terms of the above equation are of 
order 1. In terms of the original variables we get z ~ 0.85s ~ 85fim 
We can also substitute the expression for B into (J2.27|) . and we get: 



j min — Lq 1,, ,/., ,/„ , kboS 2 e z / s — kbos 2 — L Ti 



h 2 



+ -kb e z/s h = -kb e z/s s + 
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This gives us a quadratic equation in h which can be written in a simplified form by 
introducing (3 = kb e z ^ s . 

,2,o/ 1 L m i n + kbo S 2 . L m i n — Lq 

h +2(* + — ^" 7 )/>+— ^ = 

We know the numerical values of all terms in the equation, so we can estimate h. 
Conclusions: 

(1) As before, experimental data does not exist for z. 

(2) Values of h can be estimated for various values of L . (For example, for L = 0.35 
percent, h ~ 35/im.) In our estimates there is a dependence of h on L min — L which 
could be tested by running experiments with various species that have different pref- 
erences for L min , the minimum tolerable oxygen concentration. 



Second case: L falls under L min . We expect that there is no bacterial band developing 
in this scenario since the outside oxygen concentration, L , is below the minimum 
preferred concentration. There is only one region in this case, region I, from x = to 

x — z. 

In region I the equation for the bacterial density is: 

-vr' - frir + f lr l = 0, vl' + far - fal = 0; r = I =~ e~ x ' s 

B is density of bacteria at the boundary, b = Be~ x l s . The equation for oxygen is 
given by: 

L" = kb = kBe' x/s L^x) = kBs 2 e- x / s + c x x + c 2 
We use that the oxygen concentration outside the capillary is L : L(0) = L . 

L^x) = kBs 2 e~ x/s + Cl x + L - kBs 2 

To find B, c\ and z we need the following conditions: 
(1) Continuity of bacterial density at z: Be~ z l s = b 
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(2) Continuity of flux of oxygen at x = z: 

kBse- z ' s + (kBs 2 )/z(e~ z/s - 1) + L /z = 

(3) Boundary condition at x = z: Li(z) = 

By making the appropriate substitutions we arrive at: 

kb sz + kb s 2 (l - e z/s ) = L 
Again, recalling ( = z/s, a simple form of the expression is 



p-C=i-ir 

kb s 



2 



Here we must use a numerical estimate as well, since we again expect the terms to 

be the same order of magnitude. 

Conclusion: 

(1) A numerical estimate for z is possible, but, as mentioned above, it cannot be 
compared to experimental values. 



This concludes the steady state solutions. Now we can examine how the solution 
behaves on the time scale over which the band develops. 

Quasi steady state solutions 

In the previous section we provided an analytical solution for the steady state in three 
different regions: in front of, inside, and behind the band. In this section we obtain 
an analytical solution for the time interval 50 seconds to 3 minutes, during which the 
bacterial band forms, shown on Figure 12.201 We cannot obtain transient solutions 
for our system, but we know that the steady state solutions will be valid only after 
a very long time period (hours). On the other hand, the experimental observations 
are for a much shorter time scale, on the order of minutes, so we would like to find 
the quasi steady state solutions for this time period. We will show below that during 
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this time period only the bacteria in front of the band have time to move inside the 
band, so in fact, we only have two regions: one in front of the band and one inside 
the band. We check our heuristic analytical results against the experimental findings. 

b io 

9 
8 
7 
6 
5 
4 
3 
2 
1 

0.5 1 1.5 L 

Figure 2.20: Quasi steady state solution of the bacterial density. 

In order to demonstrate that only the bacteria in front of the band are involved in the 
band formation, we must estimate the time for the cells from behind the band to move 
into the band by pure diffusion in the absence of oxygen. The diffusion coefficient is 
given by: D ~ 4 - - 8 • 10 2 ^ ~ 8 • l(T 4 ^. Now we can estimate the 

° J / 0.5s 1 s s 

time of diffusion, t ~ -tt — f *l 5 ! r " ra ' > 2 , ~ 2.8 ■ 10 3 s ~ 47 min. This confirms that if 

D 8-10 *mm^/s 

we only want to look at the time interval up to 3 minutes, we only need to consider 
two regions, region I, from x = to X\ = d, which is in front of the band, and region 
II, from Xi = d to x 2 = d + h, inside the band. 
Region I: 




CHAPTER 2. 



MATHEMATICAL MODELS IN BIOLOGY 



54 



The equation for oxygen in this region is: 



dx 2 



with boundary conditions: 



L(0) = L (2.29) 
L(d) = L max (2.30) 

Using ()2.29|) we arrive at the equation for the oxygen in region I: 

L(x) = L — c\x 

We assume that there are no bacteria in this region because they have all aggregated 
to the band. 

Region II: 

The equation for oxygen in this region is: 

dx 2 



With boundary conditions: 



L(d) = L max (2.31) 
L(d + h) = (2.32) 



Using (|2.31j) the equation for oxygen in this region is 

L = L max -d(x-d) + - df 

The number of bacteria in the band is constant, b = B. 
We are trying to find four parameters, ci,d,h and B. 
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ci = - (2.33) 

a 

L max ~ Cl h + ^h 2 = (2.34) 

L'j(d) = L'jj(d) (2.35) 

L' n (d + h) = (2.36) 

Bh = b (d + h) (2.37) 

We obtained (12~33|) using IjCTljl . and (12~33J) from (|232) . In (1230) the conservation of 
flux at the left edge of the band is a redundant condition, —c\ = —c\ + (d — d). The 
expression (J2.36)) gives —C\ + ^§-h = 0. Equation (J2.37)) follows from the conservation 
of bacteria. Solving (|2.33|) and (|2.36|) for c\ and setting them equal gives L °~^ max = 
*§-h. We can solve this for d, 

D(L - L r 



d 



■'max j 



X Bh 

From (|2.34|) . again using (|2.36j) we arrive at L max — ^h 2 = which we can use to 
obtain 

, l^DL max 

From (|2~37f) we obtain B = b ^. 

Now we introduce a new variable, k = |j. Using the new notation we can rewrite as 
^ = ( ' L °~ k ^ ax ^ — kS^d+h) • ^ e i m P ose t ne condition that d, the distance in front of the 
band must be much larger than the width of the band, d » h which implies that 
d ~ -j^. This gives 



if L Q = 0.2. If L = 1, then d ~ 17. 

The expression d ~ also can be obtained in a more intuitive way. According to 
our assumption, the oxygen flux is equal to the oxygen consumed, or: D^f = xBh. 
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This implies = *§-h = kBh = kb ^h (by using the expression for the conservation 
of bacterial cells). From this we arrive at d ~ J ^ 

Our assumption on d allows us to rewrite B ~ 6 r ■ Then we can substitute this 



2.D Lfjiax r^j / 2Zy 



expression for B into the expression for h we had earlier, h __ y ,, _ y , , ; , , 
Solving this quadratic expression for h we obtain h ~ 2D ^£ X — 2 kb 7 = 2 fc^T 
For L = 0.2 this gives ft, = 0.4, and for L = 1, h = 0.2. 

Since we have an expression for both d and h in terms of the original variables, we 

can find C\ from (|2.HHjl and B from B ~ 6of 

Conclusions: 

(1) The distance of the band from the meniscus, d, is given as a function of the outside 
oxygen concentration, L : 

i- L " 



kb d 

This formula provides estimates d ~ 0.8 mm for Lq = 0.2 and d 1.7 mm for Lq = 1. 
These values are in agreement with the figures from the Zhulin experiment. (See 
Figure EHZJ) 

(2) The band width, h is given as a function of L and L max . More experimental data 
would be necessary to determine whether the dependence of h on L max is correct. 
The two values for the band width, h = 0.4 mm and h = 0.2 mm for L = 0.2 and 
Lq = 1, respectively, are also in agreement with empirical findings. 
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2.5 Conclusions 

Our main accomplishment is the development of a mathematical model for aerotaxis 
in a steep attractant gradient for fast-adapting bacterial turning frequencies. As 
discussed in Section 12.2.41 most of the literature is focused on shallow attractant 
gradients because this allows approximations and the use of the Keller-Segel-type 
chemotaxis equations ()2.1|) . As we note, the only possible approach to steep gradients 
and slow adaptation is Monte-Carlo simulations. 

The significance of our model for biologists is its prediction of the appropriate pat- 
tern formation by assuming novel receptor mechanics and fast adaptation, while the 
same patterns are impossible to obtain by models of conventional chemotaxis. This 
indicates that aerotaxis (and probably other forms of energy taxis) uses a distinct 
biochemical mechanism to achieve motility. The separation of biochemical pathways 
supports the notion [73] that energy taxis is a "flight response" which can overrule 
chemotactic behavior thus giving organisms with aerotaxis a selective advantage over 
purely chemotactic bacteria. 

An immediate benefit of our model are its estimates for the maximal and minimal 
tolerable oxygen concentrations, L max and L min , respectively. The analytical solutions 
also show dependence of the band width and the distance of the band from the 
meniscus on other parameters which makes our predictions experimentally testable. 
A possible limitation of the model is the assumption that adaptation has to be faster 
than the length of the run. A further question to investigate is the dynamic of the 
model if the run length and adaptation occur on the same time scale. The option that 
adaptation is longer than the time scale of the run is ruled out by our Monte-Carlo 
simulations. Determining whether there is an optimal characteristic adaptation time 
for given attractant gradients would also be scientifically valuable. 
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Chapter 3 

Growth cone guidance 

3.1 Introduction 

Sensory input from the environment and responses to such inputs are characteristic 
of all animal cells. Chemotaxis, the movement along effector gradients toward sources 
of attractants (e.g. nutrients or certain guidance molecules) or away from chemical 
repellents, provides an example of such a process. In fact, even bacterial cells are 
capable of chemotaxis (see Chapter |2J), although the mechanism underlying chemo- 
taxis of eukaryotic cells and of bacteria are very different. Bacterial cells, because 
of their extremely small size, are unable to detect concentration differences across 
their body. They swim through gradients, and compare concentration differences in 
time. Animal cells, on the other hand, are thought to genuinely measure concentra- 
tion differences at different parts of the cell body. Yet another difference is that the 
signal transduction pathways involved in bacterial chemotaxis of some species (e.g E. 
coli) are very well known, and accurate mathematical models of bacterial chemotaxis 
exist as well |4j. Signal transduction pathways involved in chemotaxis of animal cells 
have not been so completely uncovered, and it is possible that there are significant 
differences between the pathways as well as in the mechanisms of gradient sensing 
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between neutrophils, growth cones or slime molds (Dictyostelium discoideum). 
Understanding gradient sensing is an interesting theoretical question in its own right, 
even without considering potential biological applications of successful theoretical 
models. Many different hypotheses can be formulated regarding how concentration 
differences are compared in a cell. Dallon & Othmer ^H] discuss several distinct 
mechanisms, such as comparing receptor occupancy at different sides of the cell or 
setting up an internal gradient in response to external gradients. Theoretical models 
help to distinguish between these, and other possibilities. Many authors |IB1 EH 
also believe that the amplification of the external signal is of key importance in 
understanding chemotaxis. It is not well understood how a 2-3 % concentration 
change over its body in a noisy background is sufficient to unambiguously orient a 
cell in a gradient. While adaptation is known to be an important feature of bacterial 
chemotaxis, it is unclear whether adaptation occurs in all animal cells. 
This work focuses on chemotaxis of growth cones, partially because of the biological 
significance of the question, and partially because of exciting experimental findings 
in growth cone chemotaxis discussed in Section 13.2.11 Growth cones are the highly 
sensitive and extremely dynamic tip of axons. They are finger-like protrusions that 
are able to detect ligand gradients, and initiate the axon movement in them. Growth 
cone guidance is most important during development when axons can cover enormous 
distances to their final place in the nervous system, and also during axon regeneration 
following injuries. A successful theoretical description of growth cone chemotaxis 
could lead to very important medical applications. Growth cones are also unique 
among animal cells, because they can respond to gradients of the same chemical 
differently, depending on the receptors expressed in the cell, and also, depending on 
the internal chemical state of the cell. How the internal chemical state might influence 
growth cone guidance is explained in the experiments with Xenopus spinal neurons, 
described below. 
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The next section provides the necessary biological background, and explains experi- 
mental data relevant to the models formulated later. The second part of the Back- 
ground, Section 13.2.21 describes previous theoretical models of chemotaxis of eukary- 
otic cells. Two papers, one by Meinhardt [IB], and the other one by Levchenko & 
Iglesias are in the focus. Next, we formulate two mathematical models of chemotaxis, 
and provide analysis and numerical results of the models. The results are interpreted, 
and the we discuss the merits and shortcomings of both models. New experimental 
data on adaptation of growth cones [66 j and the limitations of the two models suggests 
some directions for further work in this area. 
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3.2 Background 

3.2.1 Biological background 

Our project aims to create a mathematical model of early events in growth cone 
response to extracellular effectors, because these events have a clear significance in 
axon guidance. During development, axons must find their appropriate targets in the 
nervous system, and they accomplish this with the help of their growth cones. Growth 
cones are able to detect spatial ligand gradients, and turn toward the source of the 
attractant. The axon body follows the guidance of growth cones in this direction. How 
growth cones are able to 'measure' gradients, even in cases where the concentration 
differences are very small, is not well understood, and is currently an important 
research topic. 

In spite of extensive research, there are many obstacles in understanding chemotaxis. 
Gradient sensing involves an enormous number of signaling pathways, and some of the 
pathways are unique to certain organisms. Many years of research has been devoted to 
studying Dictyostelium cells, for example, while growth cone signaling pathways are 
a relatively new topic of investigation. The hope is that the underlying principles are 
the same in all eukaryotic cells, however, important exceptions exist. It is generally 
accepted j^E] that chemotactic cells adapt to persistent stimulus, i.e. certain signaling 
events occur transiently shortly after the stimulus is changed. This appears necessary 
in order to explain how cells can orient in a wide range of external stimulus. However, 
it is not clear whether adaptation occurs in growth cones, as in vivo they navigate 
in attractant concentration gradients that may not span many orders of magnitude. 
Although this question is not resolved, experts in the field believe growth cones adapt 
too (Mu-ming Poo, personal communications). 

We will focus our attention on growth cone responses in netrin-1 gradients, because 
netrin-1 is the best characterized molecule which has proven to exert a tropic effect 
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in vivo. Part of the problem of gradient detection is the determination of the signal 
transduction events immediately following ligand binding, and the reconstruction 
of the signal transduction pathway. We focus on the early events, because it is 
well established that in animal cell chemotaxis sensing and motility are independent 
processes [HE], and we make the assumption that the decision to turn toward or 
away from a substance is made at an early stage of the pathway. There is also 
strong evidence suggesting that the decision is made by the time a common second 
messenger, cAMP is produced. We do not consider events happening downstream 
from this point, although we recognize that many other downstream events must 
take place for the actual turning response to be executed. 

Song and Poo [5BJ EZ] note that there are two categories of guidance cues, type I 
and type II. Members of the first group are characterized by the termination of the 
turning response when extracellular Ca 2+ is depleted, and by the co-activation of two 
pathways: the PI-3 kinase and PLC-7 pathways. Whether the turning response is 
attractive or repulsive, depends on the level of cAMP activity. In contrast, type II. 
guidance cues are independent of extracellular calcium and PI-3 kinase. The level 
of cAMP activity does not alter the turning response in guidance molecules of this 
group. The experiments of the Song, Poo, and Tessier-Lavigne labs |1H] 151] 157] 155]. 
demonstrate that turning response of growth cones toward a netrin-1 source strongly 
depends on internal cAMP levels. Several results [221 EZI indicate the dependence 
of the turning response on internal calcium dynamics as well, therefore netrin-1 is 
a type I guidance molecule. Whether the turning response to netrin-1 is attractive 
or repulsive also depends on the receptor expressed in the cell. DCC is a receptor 
that leads to attractive turning while UNC-5 leads to repulsion. However, in the 
same cell the turning response can be altered by changing the cytosolic cAMP con- 
centration. For example, in a neuron expressing DCC receptors attractive turning 
can be changed to repulsion by lowering cAMP amounts. Some of the experiments 
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investigating the calcium and cAMP dependence of the response are described in fur- 
ther detail, because we hypothesize that understanding calcium and cAMP dynamics 
leads to understanding how attractive or repulsive turning decision is made in growth 
cones. Additional information is provided about signal transduction components, in 
particular, about the interaction of the calcium and cAMP pathways. 
A number of experiments were performed whose results provide constraints for a 
model of how ligand binding leads to a turning response. The basic idea, that high 
concentrations of cytosolic cAMP correspond to attractive turning, and low concen- 
trations to repulsive turning, has been widely accepted, and is reviewed by Song and 
Poo, jHZl- This relationship is clearly demonstrated by the experiments of Ming et 
al. on Xenopus spinal neurons [48] which show that attraction to a netrin-1 gradient 
can be abolished by bath addition of Rp-cAMP, an antagonist of cAMP which causes 
lowered cAMP levels in the cell. Figure 13.11 provides information on the relevant 
parts of the signal transduction pathway. 

Song and Poo also note [HZj that for each guidance cue, there is a characteristic range 
of cAMP which determines whether attractive or repulsive turning is exhibited. When 
cAMP activity falls below the critical range, repulsive turning is observed and when 
cAMP activity is above the critical range, attractive turning is induced. The critical 
range of cAMP must be low for guidance cues that induce attraction, but further 
reduction of the cAMP converts attractive to repulsive turning. 
The justification to talk about a "critical range" or "threshold" comes from experi- 
ments done on Xenopus spinal neurons by Ming et al. [IB] • They want to determine 
whether the cAMP concentration changes turning behavior gradually, or if at a cer- 
tain concentration attractive response changes to repulsive response (and vice versa). 
They test this by administering an increasing amount of Rp-cAMP to the bath in 
addition to the same netrin-1 gradient. What they find is that there is a clear tran- 
sition from attraction to repulsion consistent with the idea of the "switch". On the 
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Figure 3.1: Part of the signal transduction pathway in growth cones. Abbreviations: AC: adenylate 
cyclase, Ca: cytosolic calcium ion, Ca 2+ ,Ca er : calcium ion in the endoplasmic reticulum, cAMP: 
cyclic AMP,DCC: Deleted in Colorectal Cancer, ER: endoplasmic reticulum, IP3: inositol 1,4,5 
triphosphate^: ligand, netrin-l,PKA: protein kinase A. References: pil5ll51 l3^1II51l4^1ISll63] 

other hand, enhancing the effectiveness of cAMP by adding Sp-cAMP to the bath (a 
chemical which enhances the effect of cAMP) has no effect on the angle of turning. 
The turning angle is statistically the same during attractive and repulsive turning 



Ming et al. also examine the effects of low netrin-1 concentrations, and they find that 
for sufficiently low concentrations there is no turning response. These experiments 
give us the following important pieces of information regarding growth cone turning 
in a netrin-1 gradient: (Note: a concentration of 5-10 at the pipette which results 
in about a thousand times the concentration at the cell leads to clear turning response 
while 0.5 and 1.5 ^ at the pipette leads to no response.) 
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• There is a certain low concentration of netrin-1 which is insufficient to elicit any 
sort of response from the growth cone. 

• Under normal circumstances, netrin-1 binding to the receptors induces turning 
toward the higher netrin-1 concentration, but lowering cytosolic cAMP levels 
can change the attractive response to repulsive turning. 

• Experiments with other molecules (rMAG) which induce repulsion, show that 
this turning response can also be reversed and changed into attractive turning 
by increasing the cytosolic cAMP levels. 

• The turning response is a 'switch-like' behavior as demonstrated by the fact that 
the angle of turning cannot increase by increasing the internal cAMP levels, and 
cannot be decreased by decreasing the cAMP levels, rather, decreasing cAMP 
levels changes attractive turning into repulsive turning abruptly. 

The abrupt, switch-like change in the turning response suggests that the internal 
cAMP levels must reach a value that we can call the 'threshold', above which level 
the growth cone is able to turn into the gradient, and below which the gradient is 
repulsive. There is no empirical evidence on how cAMP is distributed inside the cell 
during attraction or repulsion. The reversal of behavior is achieved by bath additions 
of substances which are known to lower or increase cAMP levels. If we assume that 
cAMP levels are uniform in the cell, and the attraction or repulsion is only dependent 
on whether the cAMP levels are above or below threshold, then quite clearly, we 
need another mechanism to provide the spatial information about which direction 
the growth cone must turn. On the other hand, if we assume that cAMP is elevated 
to different levels in the cell, and this is the only mechanism that drives the turning, 
we arrive at a contradiction again, because raising the internal cAMP levels by bath 
addition of Sp-cAMP should (depending on the amount of Sp-cAMP) increase the 
cAMP level everywhere in the cell above the threshold level which would terminate 
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any type of response. Based on this heuristic argument, it is unlikely that cAMP 
activation alone could be responsible for the observed behavior. Let us examine what 
is known about the Ca 2+ pathway. 

First of all, many experimentalists note that depleting the external calcium com- 
pletely, any sort of turning response is abolished [H2JIHI]- Hong et al. |S2] examine 
how different calcium levels regulate the turning response in Xenopus neurons, by 
blocking plasma membrane channels, channels of the internal storage, and by trig- 
gering calcium efflux from the internal storage by bath addition of ryanodine. They 
observe that high level internal Ca 2+ signals (resulting from normal extracellular cal- 
cium concentrations) in a netrin-1 gradient lead to attractive turning, while low level 
Ca 2+ signals (resulting from a depleted extracellular calcium level) lead to repulsive 
turning. In these experiments the authors cannot indicate whether the Ca 2+ level 
elevation is local in all growth cones, although in the larger ones they are able to 
observe a transient calcium gradient on the side facing netrin-1. 
A natural question, given the role of calcium as an intracellular messenger, is the 
spatial distribution of calcium during growth cone turning. Both Hong et al. [321 and 
Zheng [77j conduct experiments on Xenopus neurons to show that spatially restricted 
calcium elevation on one side of the growth cone is sufficient to trigger a turning 
response. Zheng is able to quantify the amount of calcium present by his technique of 
releasing caged calcium with a laser beam, and he is able to calculate the amount of 
caged calcium by controlling the frequency of the laser. He shows that the base level 
of cytosolic calcium plays an important role, because the same amount of increase 
induces attraction in a cell which has 'normal' cytosolic calcium levels, and repulsive 
turning in a cell which has lowered cytosolic calcium levels. Hong et al. manipulate 
calcium induced calcium release (CICR) by creating a ryanodine gradient across the 
growth cones. Ryanodine is known to activate the receptors of the internal calcium 
storage, the endoplasmic reticulum, ER. The calcium release from the stores is auto- 
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catalytic. Hong et al. do not quantify the amount of released calcium, only note that 
repulsive turning corresponds to shallower gradients and a lowered cytosolic calcium 
level. 

The interaction between the cAMP and calcium pathways is examined by Hong et 
al. They induce attractive turning by bath addition of ryanodine, which leads to in- 
creased calcium release, and subsequently override the attractive response by blocking 
a downstream component, PKA, of the cAMP pathway. Similarly, they can overturn 
a calcium induced repulsion by bath addition of Sp-cAMP. 

There is also a steady calcium elevation in some growth cones 10-15 minutes after 
the netrin-1 gradient is set up, which only terminates when the netrin-1 gradient is 
no longer present. This was not observed in all growth cones, and it is possible that 
this is a result of other processes (for example actin dynamics) that are unrelated to 
the initial decision making regarding the direction of turning [32J. 
We can summarize the relevant facts about the calcium pathway as follows: 

• Influx from the extracellular medium is important, because depleting the calcium 
here abolishes the turning response. 

• Local elevation of calcium mediates a response, but whether it is attractive or 
repulsive, depends on the absolute amount of calcium, not just the gradient of 
calcium across the growth cone. High levels of calcium lead to attractive turning 
in the side of the cell where calcium levels are elevated, but elevation of calcium 
to a only a low level leads to repulsive turning from the side with the elevated 
calcium. 

• cAMP can change the turning response initiated by a calcium gradient. If the 
growth cone is responding by repulsion, elevating the cAMP levels will change 
the response to attraction, and vice versa, attraction induced by high calcium 
levels on one side can turn into repulsion, if cAMP is lowered. 
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Calcium and cAMP are common second messengers and are known to interact with 
each other several ways. cAMP is produced by adenylate cyclase, and different adeny- 
late cyclase isoforms can be activated or inhibited by calcium. Adenylate cyclase in 
Xenopus spinal neurons is calcium-activated (Y. Gorbunova, personal communica- 
tion). Calcium activates three types of AC, namely AC1, AC8 and AC3, but currently 
it is not know which of these AC isoforms is found in Xenopus spinal neurons. Each 
of these isoforms responds to different calcium levels. AC3 only appears activated by 
un-physiologically high concentrations of calcium. AC8 also responds to calcium in 
concentrations 5 to 10 times higher than the concentration necessary to activate AC1, 
and is considered a "pure calcium detector" [T|J|- AC1, on the other hand responds to 
simultaneous activation by G-proteins and calcium concentrations in the physiological 
0.1-1 jj,M range and is considered "a coincidence detector". This is consistent 
with experimental data on Xenopus neurons, therefore we consider interactions of 
AC1 and calcium in our model. Also in favor of AC1 is the fact that this particular 
isoform is abundant in the central nervous system, mainly in the brain, particularly 
during development when growth cone guidance is especially relevant [To] . 
Another known effect connecting the two pathways is the phosphorylation of IP3 re- 
ceptors by protein kinase A or PKA (downstream product of cAMP) which leads 
to changes in calcium flux from the endoplasmic reticulum, ER into the cytosol. In 
some systems it increases calcium flux, in some others the phosphorylation leads to 
a decrease. We hypothesize that in our case the flux increases, and also, that phos- 
phorylation due to PKA to be proportional to adenylate cycles. This concludes all 
the biological information necessary to understand the development of our theoretical 
model of axon guidance. 
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3.2.2 Theoretical models of gradient sensing 

Theoretical descriptions of how cells can detect sometimes only a few percent change 
over their body and orient themselves toward attractants have a long history. Pre- 
vious models (Moghe & Tranquillo (50], Tranquillo & Lauffenburger [73]) focus on 
a phenomenological description of chemotactic movement. Although these models 
do include one intracellular messenger, their goal is not the accurate description of 
the signal transduction underlying gradient sensing, rather, they attempt a com- 
plete description of leukocyte chemotaxis from sensing to directional movement. Al- 
though these models made important contributions to the understanding of chemo- 
tactic movement, their approach to model chemotaxis is significantly different from 
later models which separate sensing from motility, and build on detailed informa- 
tion about signal transduction events. Since these models, separation of sensing and 
motility has emerged as an important principle in understanding chemotaxis [56J, dis- 
cussed below. Although we do not describe the Moghe & Tranquillo and Tranquillo & 
Lauffenburger models in detail, they are similar to models discussed in Section 12.2.41 
Recently, several sophisticated types of theoretical models of the signal transduction 
events underlying bacterial and leukocyte chemotaxis have been developed |U E3 El 
H31 • However, chemotaxis in eukaryotic cells is still not well characterized for a num- 
ber of reasons. One problem is the sheer number of connections and pathways that 
exist in these organisms, which makes the design of a tractable model of all signal- 
ing molecules involved in gradient detection and motility impossible with traditional 
methods. Another problem is the difficulty in deciding which experimental data is 
widely applicable to all chemotactic cells, which data may be true for Dictyostelium, 
for instance, but not growth cones. An example of such a case is the question whether 
growth cones adapt. Developing and maintaining polarity may not be important for 
growth cone chemotaxis either. Difficulty can also arise from the correct interpreta- 
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tion of the experiments. For example, Meinhardt 03] bases his model partially on 
observations about the dynamic nature of the membrane protrusions of Dictyostelium, 
however, later these protrusions were proved to be unnecessary for sensing. 
So far no theoretical model of chemotaxis in growth cones has been developed, in 
spite of the intriguing recent data on the signal transduction mechanism presented 
by Ming et al. [IB] and Song & Poo jHZ], among others. 

Parent & Devreotes jSS] established a widely accepted characterization of chemotaxis. 
Their work provides important criteria for every model of chemotactic sensing and 
movement must meet. 

• Extreme sensitivity and the ability to detect a concentration difference of as little 
as 2 % between front and back of the cell in a range of absolute concentrations; 

• Polarity: when the cell is exposed for a period of time to the same gradient the 
rear becomes less sensitive; 

• Directional sensing is not essential for movement; 

• Movement is not necessary for sensing (i.e. it is not like bacterial chemotaxis) 

• Adaptation: transient response is observed in response to uniform changes in 
the attractant concentration, while responses at the leading edge are persistent; 
eventhough uniform changes lead to transient response, immobile cells are still 
able to sense an unchanging gradient. 

Although Parent and Devreotes do not emphasize this, the amplification of the signal 
(which is necessary to explain the enormous sensitivity) is also a common goal of 
theoretical descriptions of chemotaxis. 

We review two fundamentally different theoretical models of growth cone sensing, 
one by Meinhardt |44j . and the other one by Levchenko & Iglesias [4*1] . We discuss 
how these models address the criteria set by Parent & Devreotes, and some implica- 
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tions and limitations of these models. Three other models exploring some aspect of 
chemotactic sensing are also mentioned. 

Perhaps one of the most widely quoted and most widely criticized model of recent 
years is Meinhardt's chemotaxis model jUj, which attempts to give a general frame- 
work of chemotactic sensing. His model is based on two main observations. Firstly, 
that chemotactic cells are extremely sensitive and are able to detect only a few percent 
change in the attractant concentration over the cell body, regardless of the absolute 
concentration of attractants. Secondly, that sensing is a dynamic process involving 
quickly changing protrusions of the cell membrane, called pseudopods. Even when 
no external stimulus is present, pseudopods of Dictyostelium can travel around the 
cell circumference or, in other cases, happen in synchrony on opposite sides of the the 
cell. Meinhardt's model seeks to reproduce these characteristic patterns of pseudopod 
extension. Meinhardt, echoing Parent & Devreotes, assumes that sensing and motil- 
ity are independent processes. Based on the observations, he sets four criteria for 
the model: high sensitivity; sensitivity in a wide range of attractant concentrations; 
polarization of the cell adapts to changes in the orientation of the external signal; 
intracellular pattern formation continues even in the absence of an external signal. 
Meinhardt's model describes what he calls an abstract "intrinsic pattern forming sys- 
tem" that is responsible for orienting the cell. Meinhardt first addresses the question 
of sensitivity. He proposes a Turing-like mechanism, in which a global inhibitor and a 
local activator amplify a small change in the attractant concentration. The activator 
enhances its own production as well as the production of the inhibitor. Although the 
mathematical details are not given, one can assume that this is a standard reaction- 
diffusion model, in which the inhibitor diffuses faster than the activator, therefore the 
range of inhibitor is larger than that of the activator. Such mechanisms produce a 
stable pattern which does not respond to later fluctuations in the attractant concen- 
tration. The size of the response is independent of the initial attractant gradient. The 
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model, therefore, explains how sensitivity can be independent of the absolute ligand 
concentration. However, an important shortcoming of this mechanism is that once 
the cell orients itself, it is unable to respond to new stimulus, because of the stable 
pattern of the internal signaling system. Meinhardt explores several ways a stable 
pattern can be re-set, and checks the predictions of each method against experimental 
findings. 

First, if the half life of the inhibitor is longer than that of the activator, oscillations 
occur. First, the activator level peaks, subsequently the inhibitor accumulates and 
ends the activation, and this gives one full cycle of oscillation. This model implies that 
sensing is not continuous, rather, it is possible in certain time intervals correspond- 
ing to the phase when the activator levels are low. However, there is experimental 
evidence to the contrary at least in Dictyostelium. Another way to destabilize a pat- 
tern is by reducing the range of the inhibitor to a region smaller than the whole cell 
surface. This allows more protrusions to appear, but often this orients the cell in the 
wrong direction. 

Patterns produced by reaction-diffusion systems are also adjustable when an upper 
bound is imposed on the production rate of the activator. If, in addition, the acti- 
vator is slowly diffusing, then the activated region can move around the cell surface, 
but it also becomes broad, unlike the appearance of the protrusions observed. If 
the activator is almost non-diffusible, then the protrusions have some very desirable 
characteristics. Namely, a steeper attractant gradient leads to more preferential pro- 
trusions in the cell, and without the external gradient protrusions appear randomly 
distributed over the cell surface. However, in this case the appearing peaks cannot be 
shifted. In order to explain the random appearance of pseudopods over the cell surface 
in the absence of attractants, Meinhardt keeps the assumption that the production 
of the activator has an upper bound and that the activator diffuses very slowly. This 
model still does not adjust to new attractant gradients, however. To achieve this, 
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Meinhardt assumes the existence of a second inhibitor that diffuses slowly and that 
acts on a slow time scale. The second inhibitor accumulates over time where the peak 
of the activator is located, and it destroys the activator peak. This process readjusts 
the cell, and it allows the formation of new protrusions. 

The model is analyzed exclusively numerically. The equations given by Meinhardt 
are ordinary differential equations describing the evolution of the activator and two 
inhibitors at the surface of the cell broken into n sections, a is the activator, b is 
the global inhibitor with fast diffusion and c is the local inhibitor that acts on a 
slow time scale. The constants of the equation are: s - ligand concentration; several 
different functions are used here to reflect random fluctuations or external asymmetry 
in different sections of the cell surface; b a - basic production rate of the activator; s c 
- Michaelis-Menten constant for the local inhibitor; s a - saturation constant of the 
activator; r a - decay rate of the activator; r b - decay rate of the global inhibitor; b c - 
production rate of the local inhibitor; r c - decay rate of the local inhibitor. 

dai Si(a*/b + b a ) 



dt (s c + q)(1 + s a aj) 
db 
dt 



db n 

— = n^di/n - r b b 



i=i 

dc . , v 

— = b c ai-r c a (3.1) 

These equations contain the implicit assumption that the activator, a and the local 

inhibitor, c do not diffuse while the global inhibitor, b diffuses so rapidly that it is 

given by the same function in each part of the cell membrane. The term 

Sijaj/b + bg) aj + b a b Sj 

(s c + q)(l + s a af) 6(1 + s a af) s c + c* 

shows that a«, the activator is autocatalytic, and as its level increases, it saturates. 
Both b and c inhibit the production of a. The activator decays at the rate r a . a 
promotes the production of both b and c. The production of b depends on the average 
level of the activator in the cell, while the production of the local activator only 
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depends on the local activator level, a*. The mathematical treatment of the model in 
this paper is superficial, as the author draws on his extensive experience of reaction- 
diffusion systems. 

Several aspects of the model have been attacked. It is unclear what biochemical 
mechanism the proposed pattern forming system [H] represents, and which compo- 
nents of the signal transduction mechanism play the role of the activator and the 
inhibitors. Meinhardt does suggest a molecular interpretation, but he calls it "ten- 
tative" as well. Levchenko & Iglesias note that although the existence of a second 
inhibitor is a theoretical solution, it makes the model less likely to be biologically ac- 
curate jH] . Another problem is that with the Meinhardt model persistent activation 
at a given small region of the cell membrane is impossible. Instead, activated regions 
move around the circumference of the cell. This is consistent with the dynamic fluc- 
tuations of the pseudopods that are at the basis of the Meinhardt model, however, 
these fluctuations have since been shown to be unnecessary both for the adaptability 
and persistence of signaling in Dictyostelium cells |Hj. In spite of these and other 
criticism, the model is widely known, because it does provide an appealing framework 
for the processes that results in chemotactic sensing. Specifically, the model provides 
a mechanism for amplification and the readjustment of the signaling pathways, so 
the cell can reorient itself in a changing external stimulus. Meinhardt's model is also 
attractive, because it is minimal, in the sense that it explains chemotactic sensing 
with assuming the fewest possible signaling components. This makes the model easy 
to grasp and test against experimental data. There have been few theoretical alter- 
natives proposed which can so consistently and clearly explain the most important 
features of chemotactic orientation. 

Such a theoretical alternative is offered by Levchenko & Iglesias jH]. Their interpre- 
tation of the experimental data, therefore, their modeling goals, are different from 
those of Meinhardt. They believe that chemotactic signaling pathways must be able 
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to adapt to spatially uniform increases in the external attractant concentration and 
they must also be able to signal persistently when graded stimulus is presented to 
the cell. The hypothesis that chemotactic cells adapt to uniform increases in the 
stimulus directly contradicts Meinhardt's basic assumption regarding the necessity 
of dynamic pseudopods for successful sensing. This hypothesis is based on experi- 
ments with Dictyostelium and neutrophils where certain signaling components called 
phosphoinositides are shown to adapt perfectly. There are no differences in the other 
theoretical goals. Levchenko and Iglesias also want to account for the huge internal 
amplification of the external signal, the reorientation of cells when a different stimulus 
is presented, sensitivity in a wide range of concentrations, and the independence of 
motility and sensing. 

Levchenko and Iglesias believe that the most important feature of chemotactic sensing 
is adaptation to uniform changes in the stimulus, because this is what allows cells 
to orient in a wide range of external attractant concentrations. Thus, they begin by 
building a model for adaptation. The pathway they assume is shown in figure 13.21 
Both the activator, A and the inhibitor, I are activated by the signal, S. 




R* 



Figure 3.2: Illustration of the mechanism proposed by Levchenko and Iglesias. This is the mechanism 
for perfect adaptation. The signal, S activates both the activator, A and the inhibitor, I. The output 
is the activated form of the response element, R*. 

Based on this signaling scheme, one can write down the following equations. 
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d R* 

= -k_ r IR* + k r AR 

dt 

dA 

— = -k_ a A + k' a S(A tot - A) 

^ = —k-.il + klS(I tot - I) (3.2) 
dt 

The quantities, A tot and I tot refer to the total amount of activator and inhibitor 
available. Activation of A and / depend on the signal, while the inactivation is 
constitutive. R could be found by the conservation R to t = R + R*. By assuming that 
the available substrate for S is always much larger than A and /, i.e, that A tot ^> A 
and I to t ^ I, and by non-dimensionalizing the equations, Levchenko and Iglesias 
arrive at the equations 

dv 

— = —f3ir + ail — r) 
dt 

dA 

-dT = ' {a ' s) 

— = - a (i - s) (3.3) 
The two new constants are given as follows. 

oc = k_i/k_ a 

and 

p_ {k—r/k r ){k— a /k a ) 
k—i/ ki 

It is easy to see that the steady state of the active response element, r ss is given by 
the ratio of the activator and inhibitor concentration, and it is independent of the 
signal: 

ali 



a/i + (3 

The authors show numerical simulations for the adaptation of the response element. 
Next, the authors consider amplification of the signal. This is illustrated in Figure 
13.31 The goal is to amplify the production of the output, Rl*, They note that there 
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Figure 3.3: Illustration of the mechanism proposed by Levchenko and Iglesias. This is the mechanism 
for amplification. R* promotes the production of Rl*, the signaling component which is amplified 
in this scheme. Rl* is autocatalytic, because it controls the substrate for its own precursor, Rl. 

are several ways to achieve signal amplification: for example, increasing the amount 
of enzyme or the substrate. Increasing the enzyme concentration would correspond 
to increasing the concentration of R*, while increasing the substrate concentration 
corresponds to increasing the amount of Rl. The figure shows that the authors choose 
amplification by increasing the substrate. By letting Rl* control the production of 
Rl, they create what they call "substrate-supply positive feedback". The significance 
of this mechanism is that amplification only occurs when R* is turned on. The 
implication of this is better understood when the entire signal transduction pathway 
is considered. 

Adaptation and amplification occur at different levels of this signal transduction path- 
way. Let us examine how the level of Rl* changes according to this scheme. As we 
have seen in the first scheme, R* can only be activated when a signal is present, and 
R* can adapt, or return to a base level. Only when there is a signal, S, can there 
be an amplified response, Rl*. Therefore, the autocatalytic activation is dependent 
on the existence external signal, and it cannot grow unboundedly. More importantly, 
there is no need for an inhibitor to end the signal and force the signaling pathways 
to be readjustable, like in Meinhardt's model. 

By assuming Michaelis-Menten kinetics for the production of Rl*, and assuming that 
the production of Rl is mediated by another enzyme, E, the authors give the following 
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Rl* 



Figure 3.4: Illustration of the mechanism proposed by Levchenko and Iglesias. This is the mechanism 
for perfect adaptation. This figure shows the whole signaling pathway. Signaling is possible only 
when S is present. R* adapts to persisting stimulus, and Rl* is the final, amplified output. 



equations for Rl*, Rl and E. 

dRl* 



dt 



-k- 2 Rl + 



k 2 RlR* 
k M + Rl 



dRl 

~dT 
dE 

~dt 



■ -k-xRl + (Jfei + k 2 E)Rl* 
-k_ e E + k e (E tot - E)R1* 



(3.4) 



However, so far the model only describes temporal dynamics, and does not show 
how an internal spatial gradient of might develop. If signaling takes place entirely 
locally, then adaptation means that all components of the signaling pathway return 
to the same base level everywhere, and no internal gradient develops. To answer this 
question, the authors assume that all reactions described in equations 13.21 and 13.41 
(with the appropriate initial conditions) happen in n separate compartments along 
the cell membrane and, in addition, the inhibitor is allowed to diffuse. They add the 
appropriate terms, ko(Ij+i + 1^ ~ 2J 3 -) to the equation describing the evolution of 
the inhibitor in the jth compartment. (They also assume that there is not flux at 
the two endpoints.) Now, because the greatest level of inhibitor is produced where 
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the external signal is the strongest, each compartment experiences a slightly different 
ratio a/i, so the steady state of the activated response element will be different in 
each compartment. This allows the development of a gradient in the internal signaling 
system. Numerical simulations illustrate how the model works, and verify that both 
perfect adaptation and an amplified graded response are possible. The authors also 
enter a lengthy discussion of how this theoretical scheme is mapped onto the actual 
biochemical pathways of amoebae and neutrophils. At the heart of their argument 
are G-protein activated phosoinositide kinases and phosphatases. 
The Levchenko & Iglesias model shares a lot in common with Meinhardt's model. 
They both use an activator-inhibitor system where the activator is assumed non- 
diffusive and the inhibitor diffuses quickly. However, in Levchenko & Iglesias the 
production of the inhibitor is linked to the signal rather than the activator, and 
in this model the amplification is a "substrate-supply feedback" mechanism that is 
switched on only in the presence of the external signal as well. Another significant 
difference between the models is that in Levchenko & Iglesias adaptation and ampli- 
fication happen at separate levels of signal transduction. These features solve some 
problems of the Meinhardt model. Levchenko and Iglesias offer a plausible theoretical 
framework for the understanding of chemotaxis, and the authors make some experi- 
mentally verifiable predictions regarding the nature of the inhibitor and activator. 
The two mathematical models discussed so far offer the two most comprehensive 
conceptual approaches to chemotactic sensing. Some other models must also be men- 
tioned because their results illuminate particular aspects of chemotaxis. 
Dallon and Othmer ^3] developed a mathematical model which carefully analyzes 
chemotactic signals between Dictyostelium discoideum (slime mold) cells. Their novel 
approach focuses on the signal characteristics at the boundary of the cells. The 
authors use these characteristics to make predictions regarding which, out of four 
potential mechanisms, is the most feasible to act SIS db SI gnal to initiate chemotactic 
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orientation. The four mechanisms to orient cells in a ligand gradient which had 
been proposed in the literature are as follows. Spatial sensing: the cell measures the 
concentration difference or the difference in the number of occupied receptors in the 
front and the back of the cell. Differential force mechanism: the cell adhesion to 
the substratum and to other cells depends on the level of chemoattractant. Pseudo- 
spatial mechanism: the cell extends pseudopods to convert the spatial gradient in the 
attractant into a temporal gradient. Spatio-temporal sensing: the external attractant 
gradient sets up an intracellular gradient. Adaptation is possible with this mechanism, 
by allowing the internal gradient to decay if the external gradient is unchanging. ^H] 
Dallon and Othmer focus on distinguishing between the spatial, pseudo-spatial and 
spatio-temporal mechanisms. 

When Dictyostelium cells are starved, they start secreting cAMP which acts as a 
chemotactic signal to initiate aggregation of cells. cAMP diffuses, and it is broken 
down by two chemicals: mPDE inside the cell and ePDE outside the cell . Externally, 
these are the only reactions included in the model. The signal transduction inside the 
cell is based on a previous model which postulates two pathways: an excitable one 
and an inhibitory one. These two pathways regulate the production of cAMP by the 
cell. The model assumes two cylindrical cells which are homogeneous in the vertical 
direction, so the problem can be solved in the plane. After non-dimensionalization 
and the use of the conservation for mPDE and ePDE, the model consists of the fol- 
lowing equations. Outside the cells cAMP diffuses and it is degraded by ePDE. At 
the external boundary of the cells, the outward flux of cAMP is equal to the degra- 
dation due to mPDE and secretion. Inside the cell there is another reaction-diffusion 
equation for cAMP accounting for the cAMP diffusion and its degradation due to 
mPDE. At the internal boundary, the inward flux for cAMP is equal to the basal 
production, stimulated production minus the secretion. The relevant components of 
the signal transduction mechanism are membrane-bound, and their evolution is given 
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by three ordinary differential equations. 

A first set of numerical simulations examines the early phase, when only one cell is 
signaling and the receiver cell is inactive. The difference in the cAMP concentrations 
at the front and back of the receiving cell are shown for various distances between 
the cells, and for various activity levels of mPDE. Increasing the activity of mPDE, 
which corresponds to an increase in the attractant levels everywhere, results in a 
decreased difference between the front and back cAMP concentrations. This implies 
that basing the orientation on the difference in attractant concentrations at the front 
and back of the cell is not advantageous, so spatial sensing is unlikely. In order to 
use the pseudo-spatial mechanism the cell must detect the time rate of change in the 
cAMP concentration. However, the front-to-back ratio of this is essentially a constant 
during the early phase of the signaling, and it only peaks when the rates of change in 
cAMP are negative both in the front and the back. This suggests that pseudo-spatial 
sensing would also work poorly. However, the front-to-back ratio of cAMP increases 
with increasing activity levels of mPDE, and this aspect of the signal is also not 
hampered by cAMP levels increasing everywhere around the receiving cell. Based on 
the signal characteristics of the first set of simulations, spatio-temporal sensing is the 
most feasible alternative. 

Next, the authors also include the internal signal transduction. During the initial 
time frame when cells must orient, the qualitative cAMP profiles are the same, al- 
though the peaks shift. The overall time evolution of cAMP also changes. A second 
peak of cAMP appears that corresponds to the production of cAMP by the receiv- 
ing cell. Similarly to the first set of simulations, spatial sensing still does not work 
well, because the difference between the front and back concentrations is too small. 
However, under certain assumptions regarding the effectiveness of mPDE, the pseudo- 
spatial mechanism might be useful. In general, the front-to-back difference in the time 
derivatives is too small to orient cells, but the ratio of the rates does have a peak at 
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a later time which might be a signal to initiate cell movement. The authors conclude 
that a pseudo-spatial mechanism cannot be excluded, but it may control initiation 
of movement rather than orientation. As before, the front-to-back ratio of cAMP 
concentrations gives a clear signal again, which further supports the notion of the 
spatio-temporal sensing. Although the internal cAMP gradients are weak, because 
cAMP diffuses quickly, a stable gradient is established which may be amplified by 
some other part of the signal transduction mechanism. 

The most significant argument of the article is that purely spatial mechanism is in- 
effective for organisms which must orient in a wide range of concentrations, and this 
would restrict cells to navigate only in specific concentration ranges. However, this 
may be the case for nerve growth cones. It is also important to note that measuring 
the front-to-back ratio of the attractant concentrations is an effective sensing mech- 
anism regardless of the absolute concentration. Many aspects of the problem are 
particular to chemotactic sensing in Dictyostelium: the degradation of the external 
cAMP gradient, and the fact that cells themselves can change the cAMP profile by se- 
creting cAMP. Although the particular mechanism might be very different for growth 
cones, it is likely that here also a spatio-temporal sensing based on a steady internal 
gradient is at work. 

Finally, two articles ought to be mentioned that were published very recently, in 2001. 
They reflect the renewed interest in describing chemotactic sensing. Narang et al. [S3] 
attempt to formulate a model of chemotactic sensing based on an accurate description 
of what they believe to be the most relevant part of the signal transduction pathway. 
The model addresses two questions, namely, the sensitivity of the chemotactic cells, 
and why the cellular response is only dependent on the attractant gradient while 
it is independent of the absolute concentration of the attractant. The biochemical 
mechanism the authors focus on is similar to the one examined by Levchenko & 
Iglesias. Unlike Meinhardt or Levchenko & Iglesias who create a theoretical model 



CHAPTER 3. 



MATHEMATICAL MODELS IN BIOLOGY 



83 



first, Narang et al. focus on identifying the part of the signal transduction pathway 
responsible for a polarized response in an attractant gradient. They pinpoint certain 
membrane phosphoinositides that respond to uniform attractant changes transiently 
and spatially uniformly inside the cell, while they maintain a polarized distribution 
in a gradient. 

The model has four variables: active receptors, Rio; membrane phosphoinositides, P; 
cytosolic inositides and phosophates, /; and stored phosphoinositides, P s . All four 
species are allowed to diffuse. The cell is assumed to be two dimensional and disk 
shaped, but the diffusion inside the cytosol is assumed to be fast, so radial gradients 
are ignored. This simplifies the model to one spatial variable, 8, the angle between 
the leading edge and a given point of the membrane. At the heart of the model is a 
Meinhardt-type activator-inhibitor system in which membrane phosphoinositides play 
the role of activators, and the cytosolic inositides are the inhibitors. As in Meinhardt's 
model, the activator is autocatalytic and it diffuses slower than the inhibitor. The 
receptor dynamics are assumed to follow perfect adaptation, based on the work of 
Barkai & Leibler jlj. This assumption has been disproved by recent data showing that 
perfect adaptation of receptors is not characteristic of eukaryotic cells. The stored 
phosphoinositides do not play a crucial role in the dynamics of the model. 
Numerical simulations show that the pathway responds to spatially uniform increases 
transiently. This is a result of the assumed perfect adaptation of the receptors. Sim- 
ulations in a graded stimulus demonstrate a stable and amplified response of P, the 
membrane inositides. This response is also consistent with the behavior of chemotac- 
tic cells, as it is expected based on Meinhardt's analysis of activator- inhibitor systems. 
Further numerical experiments also verify that size of the response, i.e the size of the 
membrane inositide peak, is independent of the mean external gradient concentra- 
tion. The authors of the article do obtain amplification of the signal; chemotactic 
response in a wide range of external attractant concentration; and adaptation to spa- 
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tially uniform stimulus. Narang et al., however, do not resolve the question of how 
this amplified response can be readjusted. The article proposes that a calcium surge 
could result in a destruction of the phosphoinositide peak, and the authors mention 
that they have obtained promising preliminary data supporting this hypothesis. Al- 
though it is impossible to decide based on the brief description of it in the article, 
the calcium surge may correspond to the action of the second inhibitor Meinhardt 
proposes. 

Postma and Van Haastert (HE] investigate the limitations on the localization and 
amplification of intracellular responses by analyzing the diffusion of second messenger 
molecules. During chemotactic sensing, second messengers must transmit signals 
from the cell membrane to the cytoskeleton and various locations inside the cytosol. 
The speed of the signal transmission depends on the diffusion speed of the second 
messenger, however, fast diffusion leads to the loss of spatial information. This is 
the first dilemma addressed in the paper. A related question is the amplification of 
the signal. Linear signal transduction always produces shallower second messenger 
gradients than the original stimulus, therefore a strong local amplification is needed. 
Postma and Van Haastert propose a mechanism that enhances second messenger 
gradients. 

They consider two models for second messenger production. In the first scheme 
the cell is considered to be cylindrical, and second messengers are produced at one 
end of the cylinder. The molecules are allowed to diffuse and decay. In the second 
scheme, a spherical cell is considered in a linear gradient of the external chemoat- 
tractant. Production, diffusion and degradation of the second messenger occur at the 
cell membrane. The authors find that the dispersion range of the second messenger 
is given by the expression A = a/ D m /k_i where D m is the diffusion coefficient of 
the molecule and k-\ is the rate of its degradation. This expression implies that fast 
diffusing second messengers are only able to localize if their half life is short. 
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Next, Postma and Van Haastert propose a model for signal amplification. They intro- 
duce nonlinearity by assuming that a component of the signal transduction pathway 
translocates between the cytosol and the membrane. The active receptors stimu- 
late the already membrane-bound effector molecules which begin the production of 
second messengers. Then, as the second messenger concentration increases locally, 
more effector molecules are recruited from the cytosol to the membrane, resulting 
in the amplification of the original signal. The effector translocation can be con- 
sidered a positive feedback, or local activation. By depleting the cytosol of effector 
molecules global inhibition is introduced, therefore the system is another example 
of an activator-inhibitor model. The authors do not model adaptation of the path- 
way and the readjustment of the cell. Numerical simulations demonstrate that the 
diffusion-translocation model is able to amplify an external gradient about tenfold. 
However, this amplification is smaller than experimentally observed values. The mag- 
nitude of amplification of the model also depends on the gradient in receptor activity, 
as stronger external gradients are enhanced more. In shallow attractant gradients 
the model needs to be improved, and in these situations the authors assume an 
additional mechanism: translocation of another molecule from the cytosol to the 
membrane which activates the production of second messengers. The model offers 
two important new concepts: the analysis of the dispersion characteristics of second 
messengers, and that translocation of certain components of the signal transduction 
pathway can also act as a positive feedback mechanism. 
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3.3 Mathematical models and results 

Two mathematical models are presented in this section. Each model has focused 
on a different aspect of chemotactic sensing, and both have limitations in explaining 
growth cone guidance. The goals and shortcomings of each model are discussed in 
two separate subsections. Numerical and analytical results are presented. 

3.3.1 cAMP-adenylate cyclase switch 

Based on the experimental findings, the biochemically accurate description of the 
cAMP switch seemed like the most important goal because of its key role in deter- 
mining the turning response. The main feature of the model must be the switch-like, 
"all or none" response. As the turning angle remains the same regardless of the exter- 
nal netrin-1 concentration, we concluded that the size of the cAMP response should 
be independent of the size of the external stimulus. 

Experimental evidence IHH EH] also suggests that in growth cones, similarly to 
other animals cells, we can decouple the cytoskeletal reorganization and other down- 
stream parts of the signal transduction pathway from chemotactic sensing. Turning 
occurs over the period of minutes whereas the local elevation of calcium and the in- 
crease in cAMP in response to stimulus happen much quicker, on the order of seconds. 
This allows us to consider the simplified pathway as shown in Figure 13.11 
We also believed that because growth cones only encounter effector concentrations 
over one order of magnitude, they would only need to sense a gradient within a given 
range of concentrations. This implied that adaptation to a wide range of attractant 
concentrations would not need to be considered. 

Related to size of the cAMP response is the question of how a possibly very small 
spatial gradient of an attractant (or repellent) in the extracellular medium is amplified 
into a large internal gradient. The amplification must exist, because decisive response 
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can only be expected if the internal signal is clear. This implies that even if the 
receptor occupancy on different sides of the growth cones is similar, there must be 
large enough differences in downstream parts of the signal transduction pathway 
for the growth cone to turn in the appropriate direction. We assumed that the 
amplification occurs at the level of cAMP, if it is the unambiguous biochemical signal 
for turning. There are several questions related to gradient amplification, such as the 
ability to clearly distinguish noise from signal, and physical limitations on receptor 
activation in very low and very high ligand concentrations, but we did not plan to 
address these questions. 

This model is based on the dynamics of cAMP and Ca 2+ . Namely, high levels of 
cAMP correspond to attractive turning for some guidance molecules, and for these 
substances, lowering the cAMP concentration means switching to repulsive turning; 
a gradient of cytosolic calcium leads to attraction, while the same internal calcium 
gradient at a lower overall calcium level induces repulsion. We want to give a plausible 
explanation of how cAMP concentrations are lowered and increased in a growth cone, 
and include realistic cAMP and Ca 2+ interactions. 
To summarize, we wanted the model in which 

• sensing is modeled independently of motility 

• the size of cAMP response is independent of the stimulus 

• there is large internal amplification of the external stimulus 

• there is no adaptation 

• cAMP and Ca 2+ dynamics are realistic. 

A direct approach would be including the spatial and temporal dynamics of all parts of 
the signal transduction pathway represented in Figure 13.11 This signal transduction 
pathway could be described by a system of partial differential equations with the 
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appropriate rate constants and diffusions coefficients. These constants, in general, 
are difficult to find in the literature, and such a system would have a very large 
number of unknowns. Therefore, the first step in our work is the reduction of the 
system. 

We want to find the simplest mechanism that can account for a sharp switch in 
cAMP levels. The idea of a smooth external gradient of netrin-1 inducing a sharp, 
discontinuous response in the internal cAMP concentration is very closely related 
to the idea of a smooth gradient giving rise to thresholds during development, as 
discussed by Lewis, Slack and Wolpert 

In order to model this mechanism, we made further reductions in the pathway to 
be considered, and only focused on the interaction between cytosolic calcium and 
adenylate cyclase (AC), the enzyme that produces cAMP. There is ample evidence 
|14U16ll52] of a wide range of such interactions. By focusing on these, we hypothesize 
that these are the most important nonlinear interactions contributing to the switch. 
Such a simplification is based on all other processes happening on a faster time scale. 
This assertion should be checked again if other signaling molecules would be added 
to the model. 

As discussed in the "Biological background" section, we assume that the adenylate 
cyclase isoform found in growth cones is calcium-activated, and it can simultaneously 
be activated by the receptor, via G-proteins. We also consider a positive feedback 
loop on AC by assuming that the protein produced by cAMP, called PKA enhances 
the calcium flux from the cytosolic calcium stores. This assumption is based on 
phosphorylation of the receptors on the calcium stores by PKA, also mentioned in 
the "Biological background". We assume that the increase in the calcium flux is 
proportional to the concentration of AC. This implies that there is only amplification 
between AC and PKA, and also, that production of cAMP, then production of PKA 
happens on a faster time scale then the calcium- AC or calcium-PKA interactions. We 
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formulate the model for calcium and adenylate cyclase, so it is important to comment 
on how the AC level is related to the internal cAMP concentration. A large amount of 
cAMP is produced by active adenylate cyclase, and consequently cAMP is degraded. 
No other processes regulate cAMP, therefore, the concentration of cyclic AMP is 
proportional to the concentration of active AC. The amplification of cAMP due to 
its production makes the actual switch mechanism is even more dramatic than the 
results of our model show. Our system of differential equations for calcium, denoted 
by C and the activated form of AC, denoted by A is: 




3 

dC L C 2 

-dF = ko k^TI ~ kl KfTc-2 + ^c b - c) 



+ (k f + k 3 A) ■ -° {Cer ~ C) 



dA , L C m C A 




dt k n2 + L K5 + C m C 4 

5 6 

C(0) = c b 

(3.5) 

The equation for calcium dynamics draws heavily on previous models of calcium 
dynamics j^H E3 113 ED] • The first equation describes the time evolution of cytosolic 
calcium concentration, [Ca 2+ ]. As ligand binds, there is a calcium flux from outside 
the cell which saturates with increasing ligand concentrations (1). The cytosolic 
calcium is continuously pumped into the endoplasmic reticulum (ER), as shown in 
(2). The pump is believed to transport two calcium ions per cycle, hence the second 
order form j^HEZlEl- There are a number of mechanisms that maintain the cytosolic 
calcium level near the resting value, Ca&. These include passive leak between the 
cytosol and ER, the cytosol and the extracellular medium and calcium buffering. 
These mechanisms are summarized in (3). 
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The last term, (4) describes the calcium flux from the ER into the cytosol due to the 
activation of the IP3 channels. This term is similar to the analogous term in Tang et 
al [70] who show that flux from the internal storage, the endoplasmic reticulum (ER) 
in all models based receptor-kinetics can be written in the same form. It is experi- 
mentally shown that calcium has a dual role in the dynamics of the IP3 receptors: the 
initial fast increase of calcium leads to the opening of the channels, but consequently 
calcium also contributes to the slow closing of the channel. In our model there is a 
flux purely due to the direct activation of the IP3 channels which is proportional to 
the calcium concentration difference between the ER and the cytosol, C er — C. This 
flux reaches maximal value at a calcium concentration dependent on the ER calcium 
level and k a , and is small for both low and high calcium concentrations. IP3 is taken 
to be proportional to the ligand concentration. We also include additional calcium 
flux from the ER due to the phosphorylation IP 3 channels by PKA. PKA is assumed 
to be proportional to the activated form of adenylate cycles. 

The total amount of adenylate cyclase in the cell is fixed on a short time scale. The 
inactivated part of the total given in term (7) becomes activated if simultaneously 
stimulated by netrin-1 and calcium. Activation by the ligand is assumed to have 
first-order kinetics (5). The calcium activation is mediated by calmodulin, and (8) 
is fitted to experimental data [Hj. Its fourth-order form is likely to reflect that 
four calcium ions are necessary to form Ca^CaM, the calcium-calmodulin complex 
which is responsible for the activation of AC. Activated adenylate cyclase, AC decays 
linearly with rate k$ (8). 

The parameter values are shown in Table ETT1 Wherever it was possible, we used values 
that have been established in the literature, with the exception of the total adenylate 
cyclase concentration, A t . The value for K r was obtained by fitting experimental 
data in (Hj. We tried several values, ranging over two orders of magnitude, for 
the parameters which we could not obtain from the literature. In all cases, the 
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Constant 


Value 


Value in lit. 


Reference 




7 ± 

S 


6.61 i 


EH 




1 (iM 


0.1 pbM 


E3 


h 




5^ 


EH 


K p 


0.15 fiM 


0.15 /iM 


EZI 




10 i 


2 - 

Z 2 




c b 


0.1 fiM 


0.08 /LtM 


i 




fiMs 








1 JJds 






Cer 


7 fiM 


6.3 [iM 




k a 


1 


1 






2 i 

S 








1 fiM 






Cm 


20 nM 


20 /iM 


□ 


K r 


1 /«M 




used [Tl] 


At 


20 /iM 


0.02 /iM 


□ 


h 


1 i 

s 







Table 3.1: Parameter values for the calcium-adenylate cyclase switch model 

qualitative behavior of the system remained the same. The equations 13.51 were not 
non-dimensionalized. Comments on the parameter value for A t and the reason for 
leaving the model in its dimensional form are given at the end of the section. 
It is important to mention how the model might change, if other isoforms turn out 
more important in growth cone guidance. Because AC3 responds to very high calcium 
levels, it is unlikely to play any role, however, it is possible that AC8 is present. AC8 
can also activated by G-proteins as well as somewhat higher calcium concentration 
than what is required for AC1. This would change the activation term in the differ- 
ential equation for adenylate cyclase, for example to term (5) + term (6) (instead of 
the current activation term which is (5) (6)), but based on numerical experiments it 
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is still possible to find a parameter range for which we see a similar behavior as the 
one in our current model. This implies that the main mechanism of the cAMP switch 
does not crucially depend on our current hypothesis about the AC isoform in growth 
cones. 

Time evolution of [ca] and [ac] for L=0.1 ,1,10 and 20 




Figure 3.5: Time evolution of calcium and adenylate cyclase for four different netrin-1 concentrations. 
Time units: seconds 



We solved the equations numerically using Matlab. The system is in dimensional 
form, and we used the parameter values as they are given in table 13.11 without their 
units. The numerical simulation is run for 10 seconds. Solutions for various ligand 
concentrations are shown in Figure 13.51 As some parameter values are disputable, 
the following discussion on the behavior of the model is limited to the qualitative 
behavior. Numerical values are included in the discussion of the results only in order 
to make the discussion easier, but these values are not claimed to provide realistic in- 
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formation. Increasing ligand concentrations corresponds to increases in the cytosolic 
concentration of calcium and adenylate cyclase. It is clear that changing the ligand 
concentration from 0.1 /iM to 1 /iM does not change the steady state value of adeny- 
late cyclase and Ca 2+ significantly just as changing the ligand concentration from 
10 /iM to 20 /iM does not. However, there is a significant jump between the steady 
state values as the ligand concentration increases from 1 to 10 fiM. This suggests 
that the ligand concentration is the bifurcation parameter in the differential equations 
for adenylate cyclase and Ca 2+ . 



Nullclines for L=0.1, ca-, ac — Nullclines lor L-1, ca -, ac — 




012345678 



Figure 3.6: Nullclines for values of the netrin-1 concentration, L=0.1, 1 and 10 from left to right. 
The solid curve gives the nullclinc for calcium and the dashed curve is the nullcline for adenylate 
cyclase. 
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This is verified when the nullclines of the system are plotted in Figure 13.61 For small 
values, there is only one small stable steady state. Regardless of the initial conditions, 
the cytosolic AC (and calcium) concentration will approach the same small value. 
For intermediate values, we see two stable steady states, and an unstable steady state 
separating them. If we start with low cytosolic calcium and low AC levels, the solution 
converges to the same steady state as before. However, by raising the level of L, only 
the high steady state remains. Changing the concentration of the ligand even very 
slightly changes the steady state level of AC (and calcium) drastically. This provides 
the mechanism of the cAMP switch. 

Bifurcation diagram of adenylate cyclase as a function of L 

14 | 1 , , 1 , 1 




1 2 3 4 5 6 

L 

Figure 3.7: Bifurcation diagram of adenylate cyclase as a function of L 

The bifurcations are further illustrated by Figure I3~7| showing the steady state of AC 
as a function of netrin-1 and Figure with the calcium steady state as a function of 
netrin-1. In Figure 1X71 we see that increasing ligand concentrations slowly increases 
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Bifurcation diagram of calcium as a function of L 

4.5 1 1 1 1 1 




1 2 3 4 5 6 

L 



Figure 3.8: Bifurcation diagram of calcium as a function of L 

the low steady state value of AC. This continues until netrin-1 reaches the threshold 
value of L=2.3, at which point the lower stable steady state disappears and only 
the higher steady state value remains. This results in a sudden drastic change in 
the steady state AC concentration which jumps from the lower steady state value 
of approximately 1.7 fiM to the higher steady state value of approximately 12 fiM. 
This process corresponds to following the solid line representing the small steady state 
values for AC, then making the jump through the dashed line to the upper solid line 
which represents the high AC steady state values. In fact, there is a hysteresis here, 
because now decreasing the ligand value below 2.3 (moving to the left along the upper 
solid line) will not change the steady state value of AC dramatically. The high steady 
state value decreases gradually until the ligand concentration reaches about 0.6. At 
this point we drop from the high steady state value to the low one instantaneously. 
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Similar behavior is shown for calcium. 

In terms of the chemotactic sensing, we have modeled a cAMP switch that signals 
unambiguously in certain ligand gradients. We can consider a cell divided into two 
internal compartment, both of which contain the same signal transduction pathways. 
We also note that the concentration of cAMP is proportional to the concentration of 
adenylate cyclase. In some ligand gradients that contain the threshold value of 2.3 
/J.M, a small steady state value of AC is obtained at one side of the cell where the 
ligand concentration remains below 2.3 fiM, and there is a sharp jump in the steady 
state of the cytosolic AC concentration at the part of the cell membrane where the 
ligand concentration exceeds the threshold. If the ligand concentration obtains the 
threshold value, the model does satisfy the goals as follows. The size of the cAMP 
response is nearly independent of the stimulus, as the steady state value does not 
change much with the netrin-1 concentration. This is shown in Figure 13.91 where 
ligand concentrations are changed over five orders of magnitude. The only significant 
increase in the steady state values is where the system goes through a bifurcation. 
Because of the bifurcation of the AC concentration, even a very small change in netrin- 
1 can produce a very large change in the concentration of AC. The model does not 
assume adaptation, and it is built on realistic calcium-adenylate cyclase dynamics. 
It is also clear from the above discussion that the model has serious limitations. Even 
if we continue to assume that growth cones do not adapt (although recent experiments 
by Poo suggest that they do), the model fails in ligand gradients that does not contain 
the threshold value. Clearly, this is the most important shortcoming of the model, 
because this implies that chemotactic sensing in growth cones is dependent on the 
absolute concentration of the ligand, and sensing is only able to occur under special 
circumstances. This is clearly not true. 

There are a few other problems as well. First of all, as it became clear that the model 
did not properly explain chemotactic sensing, efforts to calibrate the parameter val- 
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Steady state values of ac and ca as a function of L 
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Figure 3.9: Steady state values of calcium and adenylate cyclase as a function of the netrin-1 
concentration. Initial conditions for all simulations are given by Ca(0)=Ca;, ; AC(0)=0. 

ues and to non-dimensionalize the model were abandoned. Therefore, the numerical 
value for the bifurcation parameter, and the steady state values of Ca 2+ and AC are 
meaningless. We assume the appropriate ligand gradient necessary to orient the cell 
is presented. In this case, the cell is locked onto this direction, and it will be unable 
to re-orient itself, unless the ligand concentrations fall below the value of 0.6 muM 
again. During all the previous discussion only the temporal dynamics of calcium 
and adenylate cyclase are considered, and the spatial behavior was greatly simplified 
by the two-compartment model in which the signaling components of the two com- 
partments do not interact. Diffusion of the second messengers always occurs. Even 
considering that the fast half-life of cAMP would limit the range of diffusion, there 
must be diffusion between the two compartments, leading to a diminished difference 
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between the cAMP concentration of the two compartments. Finally, this model also 
fails to explain the experimental observation that lowering the cytosolic calcium level 
changes attractive turning to repulsive turning. 

3.3.2 Adaptation and diffusion model 

The model for the Ca 2+ -AC switch does not account adaptation, i.e. for transient 
signaling in uniform changes of the attractant concentration, and for a cell's ability 
to choose new orientation in a changed attractant gradient. It also predicts that 
gradient sensing can only occur in certain ligand concentrations. In this section we 
provide a model which corrects these problems. The model is fundamentally similar 
to the Levchenko & Iglesias model which allowed for both adaptation to uniform 
increases and persistent signaling in attractant gradients. The basic concept, based 
on the Levchenko & Iglesias model is summarized below. We include this heuristic 
explanation in order to give a general idea of the approach, and the details are made 
more concrete in this section. 

Adaptation and persistent signaling are achieved by allowing the temporal adaptation 
of some response element to any given ligand concentration. This is similar to the 
Barkai-Leibler bacterial chemotaxis model f3J where the receptors are known to adapt. 
Then, by letting other components (in Levchenko & Iglesias, the inhibitor) diffuse, 
the response element is forced to adapt to slightly different levels throughout the 
cell whenever a gradient of the "inhibitor" are produced. In these cases the process 
leads to an internal spatial gradient of the response element. Thus, uniform changes 
in the attractant concentration lead to spatially uniform, transient changes inside 
the cell, whereas an attractant gradient leads to a spatial gradient of the response 
element. The conditions under which such a mechanism can function are discussed in 
this section. Our model does not consider the amplification of the graded response, 
because we assume that once an internal spatial gradient exists, its amplification can 
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be achieved in several ways downstream from sensing, and this is briefly discussed 
in 13.3.21 Although our approach is very similar to that of Levchenko & Iglesias, we 
create a minimal model in which it is sufficient to consider two components of the 
signaling pathway, and we also offer a more detailed mathematical treatment of the 
model than the original paper. 

In summary, the goals of the model are adaptation to uniform attractant increments; 
persistent signaling in gradients; and the ability to reorient the cell when new stimulus 
is presented. A model with such features must be able to sense gradients in a wide 
range of attractant concentrations. We assume motility and sensing are independent, 
and focus on the description of sensing only. 

Perfect adaptation 

Let us consider the chemical pathway represented in Figure l3~TUl We can describe the 
pathway by the following system of equations with the appropriate initial conditions 
that we impose later. 



Figure 3.10: Hypothetical signal transduction pathway. R: stored, or "recycled" substance; A: 
"activated substance" which adapts; M: "modified substance" , part of the signaling pathway that 
mediates adaptation. 



dM 



m + \{-k a (l)M + k d A) 



dt 
dA 

~dt 



rA + X(k a (l)M-k d A) 



(3.6) 
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The variables are A(t), the activated substance and M(t), the modified substance. 
(In the original bacterial chemotaxis model, R, A and M are the number of receptors 
in the "refractory", "activated" and "modified" state, respectively.) The rate k a 
depends on the ligand concentration, I and A is a large nondimensional ratio of the 
time scales of the slow and fast reactions. In bacterial chemotaxis the fast and slow 
reactions correspond to phosphorylation and methylation, respectively. In animal cell 
chemotaxis there is no reason to assume a priori that similarly, a quick biochemical 
response is followed by slow adaptation, therefore our analysis also considers the case 
when A rs 0(1). The equations assume that the production of M does not depend 
on the concentration of R, the stored substance. This is true when the concentration 
of R is very large compared to concentration of the enzymes mediating the transition 
from R to M. Another assumption is that only the activated substance, A is able to 
become "recycled", or R, and the modified substance is not. Furthermore, we assume 
that downstream effects, such as turning, depend on the concentration of the active 
substance. 

The steady state of the system is given by 

m 

A s = - 

r 

and 

_ m(r + \k d ) k d m 
S ~ k a (l)Xr ~ k a (l)r 

so the steady state of A, the activated substance does not depend on the ligand 
concentration. M s assumes that A ^> 1. The steady state of A implies that regardless 
of the ligand concentration, 1, the concentration of the activated substance will always 
adapt, i.e. return to the same value which is intrinsic to this system, as it depends 
on two fixed rates, m and r. In the simplest case k a (l) = kl. Let us assume that the 
ligand concentration l jumps to Z x at time t — 0, and we start from the steady state of 
the system at the ligand concentration / , so A(0) = A s = — and M(0) = M s = —tt- 
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The analytical solution of the system of equations is given by (see Appendix IB. II for 
details): 

M(t) ~ M 2 + (M - M 1 )e~ T f t + (Mi - M 2 )e~ rat 

A(t) ~A S + (A 1 - A.) (e- r '* - e- r '«) (3.7) 

where we have defined 

r f ~ \(k d + kh) r s ~r— , 

+ Kli 

_m _ml + (k d /kl ) 

A s — — , A\ 



r ' r 1 + (kd/kh) ' 

mL „ , m k d 1 + (kd/kl ) m k d .„ 

M = — -f, M 1 = — -f d/ , M 2 = — -f. 3.8 
r k/ r k<! 1 + {k d /kli) r kl\ 

When the ligand concentration increases, then the concentration of A grows from 

the base line A s to A% > A s on the fast time scale Tf ~ l/ r /- Meanwhile, M, 

concentration of the modified substance drops from Mo to M\. On the fast time scale, 

the levels of modified and activated substances change, but their sum is not altered. 

On the slow time scale T s ~ l/r a , the concentration of activated substance returns to 

the base line A s from Ai, while the concentration of the modified substance changes 

from M to M 2 ^ M . Now we have a set of equations (eqn. 13. 6|) that describes the 

perfect adaptation of a system to a given ligand concentration. 

Figure 13.111 illustrates how the model works. This system of differential equations 

equations and all the following ones, are solved with the Euler method. The simulation 

is run for 800 seconds, and it shows that after a transient drop in the concentration 

of the active substance it returns to the baseline level. The time step is chosen to be 

h = 0.1. As it is clear from the analytical solution, the adaptation is slower when the 

ligand concentration is small. The parameter values in the simulations are m — 0.1, 

X = 5, k = 0.2, kd = 0.2 and r = 1. k a (l) = k ■ I in these and all the following 

simulations. The same parameter values are used for m, A, k, kd and r in all of the 

following simulations, unless otherwise stated. 
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Time evolution of A for two ligand concentrations 
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Figure 3.11: Time evolution of the active substance, A for two different ligand concentrations, 1=0.1 
and 1=1. The concentration of the active substance returns to the baseline in each case. 



Spatial models 



We are interested in spatial gradient sensing, so we must look at what happens when 
both A and M are functions of time and space. In order to gain some insight on the 
behavior of the system, let us first consider an axon consisting of two compartments. 
We want to investigate how a spatial gradient of A can develop in this model. In 
each compartment the same set of reactions happens, and in addition, there is a 
flux between the compartments. We consider the case when the ligand concentration 
jumps from l to / x in the first compartment and to Z 2 i n the second compartment. 
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The equations in the two compartment model are: 
dMi 



dt 
dA x 

~df ' 

dM 2 

dt 
dA 2 

dt 



-- m + \(-k a (l 1 )M l + k d A x ) + h(M 2 - Afi) 
-rA x + X(k a (l 1 )M 1 - kdAj + k 2 (A 2 - A x ) 
-- m + X(-k a (l 2 )M 2 + k d A 2 ) - h(M 2 - Mi) 
-rA 2 + \(k a (l 2 )M 2 - k d A 2 ) - k 2 (A 2 - A x ) 



Mi(0) = M 2 (0) = A(0) = A 2 (0) = - 

r k(l ) r 



The steady state solutions are given by 



A ls = — 

A 2s 
M u = 



M 2s 

where we have defined 



m 
r 
m 
r 

mri 
Ar 
mri 
Ar 



1 + 
1 + 



rikik 



Xr 2 k p + hk s (r 2 + \k d ) 

—rikik 
Xr 2 k p + kik s (r 2 + \k d ). 
r 2 {Xk a {l 2 ) + 2ki) + 2Xk d ki 
- r 2 [Xk p + kik s ] + Xk d kik s 
r 2 {Xk a {h) + 2ki) + 2Xk d ki 
r 2 [Xk p + kik s ] + Xk d kik s 



ri = r + Xk dl r 2 = r + 2k 2 
k = k a {h) - k a (l 2 ), k s = k a (li) + k a (l 2 ) 
k p = k a (li)k a (l 2 ) 

The calculations are shown in Appendix lB.21 If li = l 2 = l, or the ligand concentration 
on the two sides of the growth cone are the same, then 



Ais = A 2s 
Mi s = M 2s = 



and if we let 



lim Mis 

A^oo 



lim M 2s 

A— >oo 



m 
r 

mr + AA^ 
~ Xk a (l) 

m k d 



r k a (l) 
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so we have recovered the steady state solution to equations 13.61 This means that 
a spatially uniform increment in the ligand concentration still results in adaptation 
in the system, and no spatial gradient of the activated substance develops. These 
calculations are also included in Appendix IB. 21 

Now let us investigate how the flux of the modified substance, k\ and the flux of the 
activated substance , k 2 will influence the steady state of the system. It is important 
to note that the results of the following calculations and the qualitative behavior of 
the system remains the same under the assumption that A > 1. See Appendix IB. 21 
for details. 

First, assume that M is non-diffusible, so ki = 0. Then 



A ls = 

mr\ Xk a (l 2 )r 2 m(r + Xk d ) 



m 



A 2s = - 



Xr Xk p r 2 



rXk a (h 



M- 



mk d 
rk a (h) 
mk d 



2s 



rk a (l 2 ) 

The steady state of the system reveals that if the modified substance does not diffuse 
between the two compartments, then it is as if the two compartment were entirely 
separated. A and M settle into the same steady state values that they would have if 
the two compartments were not connected at all. 

We do not gain additional information from assuming that k 2 = 0, because the 
qualitative behavior of the system remains the same in this case. (See Appendix 
IB. 21 ) However, assuming that k 2 ^> 1 changes the qualitative behavior. 



lim A u 

hi— >oo 



m 

lim — 

k-2^oo r . 



1 + 



r\k\k 



lim M u 

fc2^0O 

lim M 2s 



X(r + 2k 2 )k p + hhin + 2k 2 )\ 

lim A 2t 
Xk a (l 2 ) + 2ki 



m 

r 
m 



mri 
Xr 

mri 
Xr 



lk a (h)(Xk a (l 2 ) + 2k l )-k 1 kl 

Xk a {l l ) + 2k l 
k a {l 2 ){Xk a {h) + 2k 1 ) + k 1 k 
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These calculations show that the steady state of the active substance will be indepen- 
dent of the external ligand concentrations in the limit as k 2 — > oo. Therefore A; 2 ^> 1 
leads to a diminished ability of the system to respond to stimulus. The steady state of 
the modified substance still depends on the external stimulus for A ~ 0(1). However, 
if A > 1, then M ls = and similarly, M 2s = r ^ 2 ) , so in the limit as A — ► oo, 

our system again responds to stimulus as two unconnected compartments would. 
These simple calculations above suggest that the flux of M, the modified substance 
must be larger than the flux of the activated substance, A, or k\ » k 2 . A heuristic 
explanation of this is as follows. If k 2 > ki, i.e if A were allowed to diffuse faster 
than M, then regardless of the external ligand concentration the amount of activated 
substance in the two compartments becomes the same. In this case the steady state 
of M still depends on the ligand concentration, unless the reactions between M and 
A are much faster than all other rates, i.e. if A ^> 1. If A is large, then M and 
A are allowed to exchange quickly, and the concentration of the modified substance 
will depend mainly on k a (l) and kd, therefore M la and M 2s return to the values they 
would have in case of two unconnected compartments. 

Now let us consider the case that A is non-diffusive, but M is. Let us assume that 
the ligand concentration at the first compartment is larger than in the second com- 
partment, i.e. l\ > l 2 . As we had seen from eqn. 13.61 the steady state of M is 
inversely proportional to k a (l) which we are taking to be proportional to the ligand 
concentration, i.e, the steady state of M is inversely proportional to the ligand con- 
centration. This implies that M\, the concentration of the modified substance in 
the first compartment is smaller than M 2 . The flux between the two compartments 
will increase the value of Mi (and decrease M 2 ), therefore more A\ is produced in 
compartment one than A 2 in compartment two. Because the flux, k 2 between A\ and 
A 2 is negligible compared to ki, the difference A\ — A 2 is maintained. Figure EH21 
illustrates this mechanism. 
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Figure 3.12: Illustration of the two compartment model when fci ^> k^- The first figure shows the 
two compartments without any connection. In the second figure we assume that M is allowed to 
diffuse between the compartments. This results in the creation of a spatial gradient of A. 

We are also interested in how the difference between the levels the activated substance 
in the two compartments, A\ — A 2 changes with the ligand concentration. Based on 
our previous calculations, we assume that k 2 ~ 0. Furthermore, we assume that k a (l) 
is a monotonically increasing function of I, and that k a (l) > for all /. When two 
different ligand concentrations are presented for compartments one and two, we see a 
difference in the rates k a (l\) and k a (l 2 ), thus a spatial gradient in the ligand produces 
the spatial gradient of A. We notice that \A ls — A 2s \ = if k a (l 2 ) = k a (li), and this 
is the minimum value \A ls — A 2s \ can obtain. 
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The absolute difference, \Ai s — A 2s \ is bounded below by 0, and above by — . It is 
also easy to show that for arbitrarily large differences in the ligand concentrations the 
difference between A\ and A 2 approaches (We show this by taking lim fea (i 1 )_ +00 and 
lim fca( 7 2 Because of the symmetry of the expression, lim fca (; 2 )_ >00 and lim^^o 
leads to the same result.) 

We want to know how \A ls — A 2s \ depends on the difference in ligand concentrations, 
k and on absolute size of the ligand concentrations, k s . We introduce new constants: 
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a i — ^r, °2 = kiij- + Xk d ) and a 3 = Ar. 

A ls - A 2s — f(k, k s ) — ai 



a 2 A; s + f (k 2 s - k 2 ) 

We are interested in |£, how the absolute value of the difference in steady state of A 
in our two compartment depends on the ligand difference, and J^-, how the difference 
depends on the size of the ligand concentrations. 

The difference, A ls — A 2s always increases with the increasing difference in the ligand 

concentration. (This is clear from the formula A ls — A 2s as well.) Now we look at 

how the absolute concentration level of the ligand changes A ls — A 2s - 

df -a 2 + fk s 

— ai- 



2 



dk s '[a 2 k s + ^{kl-k 2 )} 
— < if a 2 > — 

In our original notation this expression means that | A ls — A 2s \ is a decreasing function 
of the sum of the ligand concentrations if 

h(r + Xk d ) > Ar 

k a {h) + k a (l 2 ) 

or k\(r + Xkd){k a {h) + k a (l 2 )) > Ar. This relationship shows that (depending on the 
explicit form of k a (l)), there is an absolute ligand concentration at which the difference 
between the activated substance in the two compartments will be the largest. Let 
us assume that we fix the difference in the ligand concentrations k a (h) — k a (l 2 ), and 
only change their sum. Increasing the absolute concentration of the ligand beyond 
the point where 

Ar 

k a (h) + k a (l 2 ) 



fci(r + Xk d ) 

will result in decreased sensitivity in sensing, and similarly, smaller ligand concentra- 
tions also result is a loss of sensitivity. The existence of a range of ligand concen- 
trations in which sensing is optimal corresponds to experimental observations. If the 
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ligand concentrations are too low, then the noise in the receptor occupancy leads to 
errors in gradient sensing. If the ligand concentrations are too high, then receptors 
are saturated, and the cell's ability to to detect gradients is compromised again. 

Time evolution of A in a uniform ligand increment 
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Figure 3.13: Time evolution of the active substance, A in a ligand increment. 



The following numerical simulations confirm the behavior of the system. The first 
figure, 13.131 shows the response of the activated substance, A in a spatially uniform 
change of the ligand concentration. The initial condition of A is A(0) = — which is 
independent of the ligand concentration. We start with our two compartment model 
adapted to a spatially uniform ligand concentration, l = 0.1. Figure 0.131 shows the 
temporal dynamics of the system in a spatially uniform ligand step to l\ = 1. There 
is a quick drop in the concentration of A, then a slow adaptation to the steady state 
level, A ss = f. 

We see similar behavior if we start with our system adapted to a ligand gradient 
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Figure 3.14: Initial condition: A is adapted to a ligand gradient. Temporal dynamics of A when the 
two compartments are exposed to the same ligand concentration. 

(which does not reflect in the initial conditions for A), and then let the ligand con- 
centration jump to the same uniform level inside both compartments. In compartment 
one the initial ligand concentration is 0.1, and in compartment two the ligand con- 
centration is 0.5, and the ligand concentration jumps to I — 1 in both compartments. 
Although the transient levels of A are not the same in the two compartments, the 
steady state levels are. Figure 13.141 shows these simulations. The values for ki and 
&2 are 1 and 0.1, respectively. The qualitative behavior of the system remains the 
same if ki = is used. In these figures the active substance in compartment one, A\ 
is always given by A(l), and A2 by A(2). Matlab, instead of plotting the value A\ 
from 1 to 1.5 and A2 from 1.5 to 2 connects A\ and A^. These numerical simulations 
illustrate that our two compartment system responds transiently to changes in the 
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A and M in a ligand gradient, t=1 000 s 




Figure 3.15: Steady state of the two compartment system in a ligand gradient. 



ligand concentration, and it settles into a ligand-independent steady state in uniform 
ligand concentrations. 

Now we must examine the other main claim of our model, namely, that it sets up 
an internal gradient of the active substance when presented with a spatial attractant 
gradient. The numerical simulation for this is shown in Figure 13.151 This figure 
depicts the steady state of the modified and the active substance when the ligand 
concentration in the first compartment is 1% = 1, and in second compartment l 2 = 
0.5. The simulation is run for t=1000 seconds. k\ — 1 and k 2 = 0.1, and again, 
k 2 = does not change the qualitative behavior. The spatial gradient of the activated 
substance is maintained while the ligand gradient remains unchanged. This results 
in the persistent signaling of the system in ligand gradients. 

Appendix IB. 31 contains notes and comments on the analytical solution of the two 
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compartment model, as well as approximate solutions to the problem if A ^> 1. 
In a general case we can assume that both A and M are functions of time and space, 
and they are both allowed to diffuse. We obtain the following equations: 
dM „ d 2 M 



m + X(-k a (l)M + k d A) + D r 



at ox 2 

8 A 8 2 A 
^ = -rA + \(k a (l)M-k d A) + D 2 — 

A(x,0) = - M(x,0) m kd 



r r k a (l) 

dA dA 

~7{ \x=0 *7\ \x=L U 

OX OX 



dM, dM, 



x=0 <~, \x=L 



OX Ox 







We provide numerical solutions to this system on the interval (0, L) = (0, 10) in 
Figures 13.171 and 13.181 Based on our analysis of the two compartment model, we 
assume D\ 3> D 2 and A > 1. The diffusion coefficient of A, D 2 is chosen to be zero, 
but as before, D 2 C 1 would also provide qualitative similar results. The initial 
conditions are chosen to be the same as in the differential equation with A ^> 1. The 
equations are solved with a FTCS (forward time center space) method. The time 
step is chosen to be At = 0.01, and the grid size is Aa; = | = 0.11. We show three 
sets of simulations verifying that in uniform ligand concentrations there is a transient 
response, and that ligand gradients elicit a persistent graded response. 
In figure l3~To1 we examine the behavior of the system which has adapted to a ligand 
gradient, but at time t=0 we present to the cell a spatially uniform ligand concentra- 
tion, 1 = 1. The middle panel shows that the cell responds very quickly, and in one 
second the both A and M are almost uniform. (The ligand concentration is equal to 
one everywhere which is difficult to see in the figures.) Finally, the panel on the right 
shows the steady state of the system at t=1000 seconds. A returns to the baseline 
value of one, and M also becomes a constant in space. 
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Figure 3.16: The figures show the temporal dynamics of the reaction-diffusion system in a uniform 
attractant concentration. The first figure shows the initial condition, the second one the system at 
t=l second, and the third one shows the steady state of the system at t=1000 seconds. 



The following two figures show the response of a cell to an attractant gradient. In 
the first one, in Figure E?. 171 the attractant concentration, I is linear, as before, and in 
the second one, in Figure ETTHl I is quadratic. 

As before, the active substance, A is spatially uniform initially, while the modified 
substance, M is inversely proportional to the ligand concentration. The system is 
shown at t=l second on the second figure, at t=100 seconds in the third figure and 
at the steady state, in the last figure. As in the two compartment model with a 
ligand gradient, in the steady state the concentration of A is proportional to the 
ligand concentration, and the concentration of M is inversely proportional to it. This 
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Figure 3.17: The figures show the temporal dynamics of the reaction-diffusion system in a linear 
attractant gradient. The first figure shows the initial condition, the second one the system at t=10 
seconds, and the third one shows the steady state of the system at t=100 seconds, and the fourth 
one at t=1000 seconds. 

spatial profile persists, representing the persistent signaling of the system in a spatial 
gradient. 

Figure 13.181 shows dynamics in a nonlinear ligand gradient. Similarly to Figure 13. 17\ 
initial conditions of A are spatially uniform, and those of M are inversely propor- 
tional to the ligand concentration. The second subfigure shows the system after t=10 
seconds, and finally, the third subfigures shows the steady state after t=1000 seconds. 
At the steady state the modified and the activated substances are quadratic. As ex- 
pected, the highest value of M corresponds to the lowest value of the ligand, and the 
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Figure 3.18: The figures show the temporal dynamics of the reaction-diffusion system in a quadratic 
attractant gradient. The first figure shows the initial condition, the second one the system at t=10 
seconds, and the third one shows the steady state of the system at t=1000 seconds. 

highest value of A corresponds to the highest value of the ligand. Persistent signaling 
is predicted again. 

Discussion 

Our model aimed at recreating a few key features of chemotactic sensing. First, 
we wanted the cell model to respond with transient signaling to a spatially uniform 
ligand concentration and with a persistent signal in a spatial ligand gradient. The 
numerical simulations for both the two compartment model (which can be considered 
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a discretized version of the reaction diffusion equations) and the reaction diffusion 
equations show, in Figures 13. 1313.141 and 13.161 that a spatially uniform ligand con- 
centration elicits a transient response from the system, but the steady state of the 
active substance is spatially uniform. Persistent signaling is also achieved in ligand 
gradients, illustrated by Figures I3.15| 13.171 and 13.181 Both transient and persistent 
signaling can also be deduced from the steady state analysis of the two compartment 
model where the steady state of the active substance depends on the difference of 
ligand concentrations. 

Secondly, the model must allow the cell to choose new orientation in a changed 
attractant gradient. This is clearly the both the active and the modified 

substance depend on the ligand concentration through k a (l), the rate of production 
of A from M. 

Finally, this model does not have the limitation of the calcium-cAMP model, as it 
allows sensing in all ligand concentrations. The adaptation-diffusion model predicts 
that there is an optimal range of ligand concentrations when k a (li) + k a (l 2 ) = fcl r^+XO • 
A fixed ligand difference will result in the largest signal when this condition is satisfied. 
The model also predicts that increasing the ligand difference results in an increase in 
the signal. 

We can also improve the model by showing the simple modifications can result in 
the amplification of the signal. In all previous simulations parameters were chosen 
so that the ligand concentration is amplified moderately by the concentration of the 
active substance, however, this need not be the case for our current system with other 
parameter values. It is also important to note that huge amplification of the external 
signal is characteristic of chemotaxis, and such amplification is impossible to produce 
with our current model. As in the Levchenko Sz Iglesias model, we assume that am- 
plification takes place downstream from A, because this assumption guarantees that 
when the external signal, I is terminated, the internal signal stops as well. The na- 
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ture of this amplification could be similar to what Levchenko & Iglesias has assumed, 
shown in Figure 1331 where A would play the role of the activated response element, 



Another possibility for amplification is based on Goldbeter & Koshland [22] who 
show that amplification can result not only from nonlinear interactions, but also 
form covalent modifications under certain conditions. Figure 13. 191 shows a reaction in 
which A is the enzyme in the production of an active response element, Rl*. Based 
on Goldbeter & Koshland, there can be a drastic amplification in the amount of Rl* 
produced when the enzymes A and E are saturated, so the total amount of free A and 
E are negligible when compared with their concentrations in the complexes produced 
with Rl and Rl*. This is one possible mechanism for amplification of A. 



Figure 3.19: A as an enzyme in the production of an activated response element, Rl*. 

It is clear both from the analytical solution to the perfect adaptation model 13.71 and 
the analysis of the steady state solution of the two compartment model that the size 
of A is very important. In the original equations 13.61 A is the non-dimensional ratio 
of the fast and slow time scales which means that A determines how fast the initial 
transient response is in comparison to the adaptation. In animal cell chemotaxis 
the time scales for a transient response and for adaptation to a signal are not well 
known. We have run the numerical simulations for the two-compartment model and 
the react ion- diffusion model for several values of A. The only observation that can be 
made based on these simulations is that for small values of A the initial transient of 
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A is longer in comparison with the time it takes for A to return to its baseline value. 
Changing A by three orders of magnitude did not appear to alter the qualitative 
behavior of the system. Further analytical treatment of the problem is necessary to 
understand the role of and the constraints on the size of A. 

Finally, we must discuss the limitations of the adaptation-diffusion model. So far we 
have not linked the activated or the modified substance of the model to any particular 
components of the biochemical pathway, because there are no clear candidates for 
such connections. It is possible that one might find chemical species that adapt 
to a certain ligand concentration, but based on the current experimental data one 
cannot say with certainty what it might be. Related to this, our original aim was to 
explain the results of the calcium and cAMP experiments on Xenopus neurons in the 
framework of chemotactic sensing. We have not addressed this, although in Appendix 
IB .41 we propose an extension of the current model that accounts for how changing the 
absolute calcium concentration changes turning behavior in growth cones. 
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3.4 Conclusions and further direction 

We presented two mathematical models that attempted to describe chemotactic move- 
ment of growth cones in response to a netrin-1 gradient. In our first model we assumed 
that the signal transduction pathways would not be able to adapt to a constant signal. 
In addition, we aimed at basing our model on experimentally observed phenomena, in 
particular, the way turning response is determined by cytosolic cAMP and Ca 2+ con- 
centrations. This model described a cAMP switch which is very sensitive at a certain 
threshold concentration of netrin-1. Mathematically, the threshold value is a bifur- 
cation parameter. The system goes through a bifurcation as the two stable steady 
states separated by an unstable steady state change to one stable steady state. This 
model did not successfully explain how gradient sensing is possible in concentration 
ranges that do not include the threshold concentration. This is a severe limitation 
of the model, because, a growth cones moving past the threshold level concentration 
would permanently lock onto the same direction, even if the ligand gradient changed. 
In addition, it is possible that the assumption of no adaptation is also flawed. 
The second model attempts to satisfy most of the criteria set for chemotactic sensing, 
except internal signal amplification. This model demonstrably explains how a graded 
internal response develops in ligand gradients, and how a uniform increase in the 
ligand concentration leads to a transient internal response. Although the model is 
theoretically more sound than the first one, it is equally limited to that, because its 
lack of connection to experimentally observable signaling pathways. The variables of 
the model may represent particular chemicals, or they might represent many compo- 
nents of the signal transduction pathway which acts as a unit on a time scale faster 
than what is considered in the model. 

Much of the model builds on the understanding how perfect adaptation occurs in 
bacterial chemotaxis. However, there are important differences between bacterial 



CHAPTER 3. 



MATHEMATICAL MODELS IN BIOLOGY 



119 



and growth cone chemotaxis which might influence the mathematical analysis of the 
problem. One such example is the known separation of time scales between phospho- 
rylation and methylation in bacterial chemotaxis which allows simplifications in the 
mathematical treatment of the problem. The same simplifications may be incorrect 
in the description of growth cone chemotaxis. 

Clearly, this topic is open for further theoretical and experimental research. Experi- 
mental observations of Song & Poo [US] on adaptation might be a promising starting 
point for further work. Song & Poo believe that growth cones periodically lose and 
regain their sensitivity to gradients. Such adaptation (which is different from the 
adaptation assumed in our second model) appears similar to an idea proposed by 
Meinhardt He states that in a react ion- diffusion system, if the half-life of the 

inhibitor is shorter than that of the activator, then oscillations occur. The accumu- 
lation of the activator is overtaken by the inhibitor, and for the period of time that 
only the inhibitor is present, the cell is unable to respond to new gradients. However, 
in order to pursue a mathematical model based on this idea, more biological data is 
necessary to formulate a hypothesis. 
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Chapter 4 

Endothelial cell deformation 

4.1 Introduction 

The present chapter concerns a topic very distinct from gradient sensing, the subject 
of the previous sections. In the model developed in this section, we investigate the 
mechanical effects of blood flow on endothelial cells. As opposed to previous models 
in which we aimed to understand biochemical signal transduction pathways, here we 
investigate a mechano-transduction pathway, i.e. a pathway which transmits physical 
signals such as forces and deformations. 

Endothelial cells form a monolayer inside blood vessels, acting as a boundary between 
the blood flow and the vessel walls. The endothelial layer is exposed to various me- 
chanical stresses, such as pressure, circumferential stretch, and tangential shearing 
forces due to blood flow. Endothelial cells respond in a wide variety of ways to these 
forces, and their response is thought to protect the arterial system from potential 
damages. Evidence supporting this hypothesis is offered by experiments which have 
demonstrated that the development of certain vascular diseases, such as atheroscle- 
rosis, coincides with the failure of the proper responses of the endothelium to flow. 
There is an array of events that takes place in endothelial cells exposed to flow, some 
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of which are immediate upon the application of shear stresses, and some, which occur 
on the time scale of many hours. Among the fast responses are activations of flow- 
sensitive ion channels and activation of G-proteins. Changes in gene expression are 
also observed, although on these events happen on a slower time scale, and finally, 
cells go through morphological changes after they have been exposed to shear stress 
approximately a day. In steady and in pulsatile flow (a superposition of oscillatory 
and steady flow) endothelial cells tend to elongate and align with the direction of the 
flow. Such enormous changes require cells to extensively rearrange their cytoskeletal 
structure. Although the cytoskeletal reorganization in response to flow induced shear 
stress has been studied, it is still an open area of research. 

The morphological changes in the cytoskeleton are well documented ( 3J original 
source: [2H1 155]). but it is unknown how the cells are able to sense shear stress, 
and once it is detected, how shear stress is transmitted from the cell membrane to 
the cytoskeleton. There have been quite a lot of previous theoretical and numerical 
investigations of the cytoskeletal changes produced by shearing forces. Our aim is 
to incorporate the effects of flow on the cytoskeleton, but in addition, we want to 
examine the viscoelastic behavior of other structures, such as the nucleus, cell-cell 
adhesions, and focal adhesion sites. However, as the cytoskeleton is thought to be the 
main force-bearing structure of endothelial cells, much of the deformation and me- 
chanical response must come from it. For this reason, we summarize some previous 
results regarding the reorganization of cytoskeleton without attempting to provide an 
exhaustive review of theoretical work in this 

Theoretical models focusing on how shear stress is transmitted often make one of 
two assumptions. They either propose that stress is transmitted by producing de- 
formations in filaments, or, that there exists an internal mechanical tension inside 
cells independently of the external shearing forces, and that instead of significant 
deformations, there is simply a rotation or change in spacing in order to respond to 
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stress. The first assumption is consistent with a model by Satcher & Dewey, called 
an "open-cell foam" model, the second assumption is used by Wang & Ingber ([69j, 
original source: [75]). and by Stamenovic et al. [UH] in so called "tensegrity" models 
as well as "cable net" models. 

Satcher & Dewey [02] investigate how shearing forces distort the polymers of the 
cytoskeleton: F-actin, intermediate filaments and microtubules. The thesis of their 
work is that F-actin stabilizes cells by decreasing deformability. The role of stress 
fibers, which are microfilaments connected into bundles, is also considered. Satcher 
& Dewey model the F-actin filaments as open lattices. (The structure of the actin 
filaments is similar to other material, such as glass foams, etc, hence the name "open- 
cell foam model".) The advantage of such a model formulation is that the Young's 
modulus, the measure of the material's ability to resist distortion, can be expressed 
in terms of filament properties, instead of having to find the density and moment of 
inertia for the entire F-actin cytoskeleton. With the open-lattice model properties of 
the cytoskeleton (shear modulus, which is the coefficient of the rigidity of a material, 
and the modulus of elastic deformation) are computed. The obtained values are the 
same order of magnitude as experimental data. The article also shows that although 
stress fibers increase the rigidity of the actin network, their elastic deformation is too 
small to effect the network, therefore their role is unclear based on the model. 
Stamenovic & Coughlin (HHl compare predictions of three different types of models 
of the cytoskeleton. The first type is the "open-cell foam" model, the second is the 
"cable net" model, and the third is the "tensegrity" model. The latter two share the 
assumption that pre-existing tensions are present in the cytoskeleton even in absence 
of external stress. In the cable net model the actin filaments are represented by elastic 
cables that are pulled tight by various forces generated inside the cell. The tensegrity 
model is similar: it assumes cables connected to rigid beams or "struts", and the net- 
work of struts and cables represents the cytoskeleton. The basis of comparison of the 
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three models is the model's estimate of Young's modulus for the cytoskeleton. This 
question if further complicated by the large range (10° — 10 5 Pa) given for the Young's 
modulus attained by different experimental techniques, such as magnetic bead mi- 
crorheometry, magnetic twisting cytometrymicropipette aspiration and atomic force 
microscopy. The paper discusses possible ways experimental procedures might bias 
the obtained values. 

The Young's modulus predicted by the open-cell foam model were much higher 
(10 3 — 10 4 Pa) than the Young's modulus predicted by the cable net and tensegrity 
models (10 — 10 2 Pa). The article concludes that different models may be appropriate 
for modeling cytoskeletal changes under different conditions. For example, the cable 
net and tensegrity models may be applicable in low stress whereas under large stress 
the cytoskeletal filaments bend, and the open-cell foam model provides a better de- 
scription. The models compared by Stamenovic & Coughlin describe only the elastic 
properties of the cytoskeleton. However, Satcher & Dewey in [62J note that the open- 
cell foam model lends itself easily to the description of viscoelastic properties, if the 
open space between the lattices is assumed to be filled with a liquid. 
Now we return to the larger goal of understanding how shear stress sensing and 
transduction occur in endothelial cells. One hypothesis is that shear stress deforms 
flow-sensitive parts of the membrane, such as certain ion-channels or receptors. The 
deformation of these structures could then be immediately transmitted to cytoskeletal 
elements connected to them. Such a mechano-transduction pathway may complement 
other, biochemical signaling pathways. A realistic model of the mechanical signaling 
pathway would allow quantitative tests of whether deformations and stresses gener- 
ated in the cell would provide a sufficient signal. A previous model by Barakat |3] 
shows that flow sensors, modeled by a viscoelastic body, respond differently in oscil- 
latory, pulsatile and steady flows, and the differences in the response could provide 
the necessary signal for downstream components of the pathway. 
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Our current work extends this model, and we represent other parts of the endothelial 
cell, namely actin filaments, the nucleus, and transmembrane proteins as viscoelastic 
Kelvin bodies as well. The final goal of this line of investigation is the development 
of a complex network of viscoelastic bodies where each body is described by exper- 
imentally obtained parameters. This dissertation is only concerned with two small 
networks, one consisting four, and the other one of seven viscoelastic bodies. Our 
work derives the equations for single bodies in series and in parallel, however, further 
work is necessary to describe more complicated connections (for example, n bodies 
connected in series which are connected to in parallel to n bodies in series again). 
Numerical simulations of the model are generated to test the model's dependence on 
parameter values, and numerical solutions of the deformation due to shear stress of 
the two simple model networks are also given. We discuss the implications of our 
results regarding endothelial cell behavior. 
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4.2 Mathematical Model 

In this section, first we derive the equations describing the deformations of coupled 
Kelvin bodies and discuss solutions to the equations. Next, we describe the networks 
we use to model endothelial cells, and finally we discuss how the parameter values 
were obtained for parts of the network. 

4.2.1 Kelvin bodies in series 

Our goal is to develop a mathematical framework to describe the deformations of the 
cell surface and intracellular structure within endothelial cells with respect to steady 
and oscillatory flow. Kelvin bodies are the most general models for viscoelastic ma- 
terials, and they have frequently been used to model how the deformation of cell 
tissues depends on the forcing [3], |20j. In addition, experimental data is also avail- 
able which describes the viscoelastic properties of the cell nucleus [24], cytoskeletal 
structures jBI] and transmembrane proteins |6] in terms of the parameters of a Kelvin 
body. This makes the Kelvin body a very effective tool for theoretical modeling of 
endothelial cell deformations. 

First, we give the equation relating the deformation and the force exerted by the 
flow for one Kelvin body as derived by Fung [271]. (The derivation is very similar to 
the case where two Kelvin bodies are connected in parallel which is shown in detail 
below.) 

The deformation u(t) of one Kelvin body as a function of a given forcing, F(t), is 
obtained by solving a first order linear equation (Fung, 20j): 



F H F = fi 01 u +771(1 H )u 



(4.1) 
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Figure 4.1: Diagram to illustrate one Kelvin body. The parameters to characterize the body are 
as follows. Dashpot viscosity: r)±, spring constant in upper branch: fin, spring constant in lower 
branch: fio\. 



The solution in steady flow, F = F is (Barakat, 3J): 



u(t) 
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(4.3) 



T e = • 

Mil 

T e is the relaxation time for a Kelvin body under constant strain (i.e. the stretch 
per unit length), and r a represents the relaxation time for constant stress (i.e. the 
force per unit area). The deformation of the two springs in the Kelvin body is 
instantaneous, as it is seen from the initial condition for the deformation. 14 .21 while the 
dashpot slowly creeps to the steady state of its deformation. Examining the expression 
for the deformation in steady flow (Eqn. 14. 3|) . it is clear that large coefficients of 
viscosity lead to longer relaxation times under constant stress while large spring 
coefficients result in decreasing relaxation time. The deformation in oscillatory and 
pulsatile flow is also given by Barakat, 

When Kelvin bodies are coupled in series, the deformation for each body can be found 
individually, so in essence it is exactly like the one-body case. The force acting on 
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Figure 4.2: Diagram to illustrate n Kelvin bodies in scries. The parameters to characterize the ith 
body are as follows. Dashpot viscosity: r]i, spring constant in upper branch: jiu, spring constant in 
lower branch: /iQj- 

each Kelvin body is the same, so the overall deformation can be calculated as the 
sum of deformations, Ui for i=l,...,n where each deformation is the function of the 
same forcing function, F: 

F + —F — tiQiU + r]i(l + —)v,i 

[I'll Uli 

M Q ) = T 

A*0i + A*li 
ii 

1=1 

4.2.2 Kelvin bodies in parallel 

We want to find the deformation for Kelvin bodies in more complicated networks. 
We begin by taking two Kelvin bodies coupled in parallel. Each of the Kelvin bodies 
is described by three parameters: the viscosity of the dashpot, and the two spring 
constants, as Figure IP1 shows below. 

In both the upper and the lower Kelvin body there are some relationships that must 
hold, namely, the total deformation, u must be a sum of the deformations of the 
dashpot and the spring in the upper branch, and this deformation is the same as the 
deformation of the spring in the lower branch. We must also note that the deformation 
of the upper body and the deformation of the lower body must be identical. These 
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Figure 4.3: Diagram to illustrate two Kelvin bodies in parallel. The parameters to characterize 
the bodies are as follows. Upper Kelvin body: dashpot viscosity: r/n, spring constant in upper 
branch: fj,u, spring constant in lower branch: ^ i- Lower Kelvin body: dashpot viscosity: 7712, 
spring constant in upper branch: [112, spring constant in lower branch: fi 02 . 

relationships give us the following equations. 



Ul2 + U'l2 



U 



u 



Another observation is that the total force, F of the two bodies splits into the force in 
the upper body, F 1 and the force in the lower body, F 2 . The total force in the upper 
body is also given as a sum of the force in the upper branch and the lower branch, 
and similarly for the lower body. This is described by the equations: 

F 1 + F 2 = F 



Fn\ + F] 



11 



Fn2 + F] 



12 



Now let us consider the upper body only. The same force in the upper branch, Fu 
is transmitted from the dashpot to the spring. The force acting on the spring is 
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proportional to the deformation it produces, and the force acting on the dashpot is 
proportional to the velocity of the dashpot. Using the variables of the diagram we 
can write this as: 

Fu = rinUn = Hx\u' n . 

Here tin denotes the derivative of un, the deformation of the dashpot, and u' n is the 
deformation of the spring. The force in the lower branch of this body acts entirely on 
the dashpot, and here the deformation is going to be the sum of the deformation due 
to the dashpot plus the deformation due to the spring in the upper branch, therefore 

F i = AioiOii + u' n ). 

Similarly, we can write down the corresponding equations for the lower body as well. 

Fl2 = ?7l2«12 = H12U' 12 
F 2 = Ai02(«12 + U' 12 ) 

Using the nine equations above, one can derive two differential equations that describe 
the the deformation of the two coupled Kelvin bodies as a function of time and the 
force acting on the bodies, with the appropriate boundary conditions. The derivation 
is similar to Fung's [20]- I n the equations below we assume that F(t) is given, therefore 
we can also find F(t). For the particular forcing functions we are interested in, the 
derivative always exists, and it is continuous. 

Fi H Fi = /i 01 u + 7711(1 H )u 

Mil /in 

u(0) = -^L (4.4) 
A*oi + Mn 



F 2 + ^^-F 2 
M12 



1 /-, . /i02x . 

/i02« + ?7i2(l H )u 

/il2 



u(0) 



F 2 (0) 

/iQ2 + /il2 



(4.5) 
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Now we can use the fact that Fi and F 2 sum to F to substitute a(t)F(t) = Fi(t) and 
1 — a(t)F(t) = F 2 (t). a(t) is the coefficient of force splitting, and it is an unknown 
function of time. 



aF + — (aF) = /a im + r] n (l + 

/ill /ill 



u(0) 



a(0)F(0) 
/ioi + /ill 



(4.6) 



;i _ a ) F + J a ) F ) = f, Q2U + ^(i + t^)u 

/il2 A*12 



«(0) 



;i-a(0))F(0) 



(4.7) 



/i02 + /il2 

Now we have a system of two differential equations with initial conditions and two 
unknown functions, u(t) and F(t)a(t), and we would like to solve for them. Solving 
the equations for u(t) and a(t) is not practical, for two reasons. For particular choices 
of the flow, for example, for oscillatory flow F(t) = Fq cos(o>t), the matrix A, defined 
below, would depend on the forcing function, and A^ 1 would become singular peri- 
odically when cos(u;t) = 0. Also, matrices A and D (also defined below) would both 
depend on time, and this could considerably slow down the computations. In order 
to avoid these problems, we compute u(t) and F(t)a(t). F(t) is a know function of 
time, so it is always possible to find a(t), if necessary. Now rearranging equations 14.61 
and 14.71 gives us the following. 

/±1V Mil 1 Mil 



7/12(1 



M02 ' 
M12 - 



r)i2 
A* 12 





ii 




-/ioi 1 




u 








-/i02 -1 







+ 



M12 



«(0) 




. o(0)F(0) _ 





F(0) 



^01+^11+^02+^12 

F(0)( m + M n) 
At01+Mll+A*02+A*12 



(4- 
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We can let 



.4 




V12 
A»12 



Vii 
Mil 



1 



D 



1 







c = 



F + 



mi p 

A*12 



U = 



U(t) 

a(t)F(t) 



Now we can write the system of equations as 




= Du + c. 



In order to compute the solution for w, we can express the system of differential 
equations in the form 



This is the equation with the appropriate initial conditions that describes the dynam- 
ics of two Kelvin bodies coupled in parallel. 

Now we can look at a generalization of the two body problem to deriving the dif- 
ferential equations governing n Kelvin bodies coupled in parallel. Let us start again 
with the diagram of the bodies with their appropriate parameters, shown in Figure 
14.41 As before, the sum of the forces in the branches has to be the total force, F. 




(4.9) 



u(0) = u 



(4.10) 



71 

i=l 



CHAPTER 4. 



MATHEMATICAL MODELS IN BIOLOGY 



132 




Figure 4.4: Diagram to illustrate n Kelvin bodies in parallel. The parameters to characterize the 
bodies are as follows, ith body (for i=l,...,n): dashpot viscosity: rju, spring constant in upper 
branch: fin, spring constant in lower branch: fi i- 



In general, we do not know how the forces split into these branches because this 
depends on the particular parameter values of the Kelvin bodies. Therefore we can 
call the force splitting coefficient for the ith branch to be Oj(i), and the force in this 
branch <ii(i)F(i). Using the above relationship we get the following: 

n-l 
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The equation for the the ith body for i — 1, n — 1 is: 



diF H -(diF) = fi iU + r]u(l H -)u 



u{0) = g.(0)F(Q) (4 n) 



And for the nth body we have a similar expression: 



71— 1 , 71—1 

(1 - V a,)F + ^ • 1((1 - V a t )F) = ^ n u + %B (1 + 

~ Uln dt yUl n 

w(0) = - ^ t=1 y !! y ! (4.12) 

Now we must rearrange the equations so we are solving for u(t) and cii(t)F(t) again. 
When we rearrange differential equations we get for i = 1, {n — 1): 

77ij(l + —)u - —(diF) = -fi 0i u + ciiF 
fJ'U [i-n 

for n: 

71—1 71—1 

ifr„(l + — + — E = - V a,F + F + 

/"in A*i» /-tin 

We also need to find the appropriate expression for the initial conditions. We have, 
first for i = 1, (n — 1): 

u{0){fi 0i + Hii) = d i {0)F{0) (4.13) 

and for i=n: 

u(0) = M>5^) ( 4. 14 ) 
A*0n + A*i ft 

w(0), the initial deformation of all of the bodies is the same, so we get n equations 
and n unknowns: w(0) and cij(O) for i=l,...,n-l. The force splitting coefficient of the 
nth body is already determined from this to be a n (0) = 1 — YllZi a i{Q) ■ We must 
rearrange the equations I4. 131 and 14. 141 to solve for u(0) and Oj(0). 

71-1 

u(0)(/i o „ + //in) = F(0) - E Oi(0)F(0) 

7=1 
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n-l 



«(0) (/ion + urn) = F{0) - «(0) (A*0i + 

1=1 

n 
i=l 

Therefore the initial conditions are: 

«(0) 

fli(0) 



Er=i(A*oi + A*ii) 
^(0)(a*w + A«i») 



Er=i(^w + ^ii) 
Now we can look at the equations in matrix form: 



A = 



.mi o ••• 

»ii 



Ml 2 



"In 
"In 



.211 

Mil 








»7l(n-l) 
Ml(n-l) 

Vln 

Mln 



D = 



-A*oi 

"A*02 



1 ••• 
1 



— A*oi 

-/iO(n-l) 
— A*0n " ' 
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Ml 2 

Just like before, we have the differential equation for u(t), annxl vector whose entries 
are u(t) and a>i(t)F(t) for i — 1, (n — 1) with the appropriate initial conditions: 

M{t) A^Du+A^c 
u(0) = u 



dt 



This is the same linear equation as 14.91 and 14.101 with the matrix A l D and vector 
A~ x c defined appropriately, and its solution is given by 



u 



{u + D- l c)e M - A- X c 



(4.15) 



where A is a diagonal matrix whose eigenvalues are the same as the eigenvalues of 
the matrix M = A~ 1 D. (Obtaining the solution to equations 14.91 and 14. 101 is similar 
to the derivation shown in Appendix IB. 31 ) Let us first discuss how to find D~ x c and 
v4 _1 c, then turn to finding A. If D~ 1 c = y, then c = Dy, in other words, 













F -f 2i« p 

Mln 





1 






1 



-/^(n-l) 
—f^On — 1 






1 

-1 





yi 








Vn-l 







In the ith row we have —fioiVi + Vi+i — 0, and this implies that for i=l,...,n-l 

Vi+i = Vim- 
la the nth row we get —fio n yi — Yll=2 Hi = F + ^-F. By substituting the expression 



for yi + \ in terms of y\ into this 
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n 



yi(-/^on - Y^oi-i) = ^ + — ^F- 

7^ 



_P _|_ Hl2L_P 



En 
i=2 /^Oi-l 

This gives an explicit formula for yi and based on our expression for in terms of y\ 
we can find all the other components of y. 

We can use a similar argument to find x = A~ x c. If we let hi = rfu(l + ^) and 



di = —— then we get that 



_P _|_ UllLp 



h n d n Yli=i dl 

_ hjXi 
di 

Now let us return to the matrix of eigenvalues, A. We are looking for the diagonal 
matrix whose entries Ai, A n satisfy the equation A~~ 1 Dx = \x for i=l,...,n and for 
i^O. This is equivalent to (D — XiA)x = which is called the generalized eigenvalue 
problem. It is easy to see that the matrix (D — \A) has the same very nice and sparse 
structure that both A nd D have, but solving the generalized eigenvalue problem leads 
to having to find the roots of an nth degree polynomial. In the most general case this 
can only be solved numerically, but some special cases of the problem can be solved 
analytically (for example, if all n Kelvin bodies have the same parameter values). 
Other methods can also be applied to solve the original system of linear differential 
equations, for example Laplace transform methods, but they lead to the same problem 
of having to find roots of an n-degree polynomial. 

In this dissertation we only use networks with two Kelvin bodies in parallel, and in 
these cases the deformation can be found analytically as well as numerically, because 
it only requires finding solutions to a quadratic equation. Extensions of the model 
to a large number of bodies in parallel would have to be done numerically, and even 
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analytical solutions can only be found by numerically computing the roots of an 
nth-degree polynomial. 

4.2.3 Model networks 

The aim of our model is to gain further understanding of how flow over the surface of 
endothelial cells leads to regulation of the cell shape. It is known that the cytoskele- 
tal structure is re-organized in a period of approximately a day, but there is no clear 
evidence to what extent biochemical events, and to what extent purely mechanical 
processes contribute to this. We want to examine how the flow-induced shear stress 
which deforms flow sensors is transmitted through the cytoskeleton to the nucleus and 
to other transmembrane proteins such as ion channels and attachments to the sub- 
strate. We want to use networks of coupled Kelvin bodies to model endothelial cells, 
and we want to investigate how deformations of the individual parts will contribute 
to the overall deformation of the cell. 

The cytoskeleton consists of three types of polymers: actin filaments, microtubules 
and intermediate filaments each of which deform if shear stress is applied to the cell 
|62j . Qualitatively, actin filaments rupture at relatively low strain, but actin can be 
rapidly recycled and filaments re-formed as it is required for cell motility, among its 
other functions. Below a critical strain actin networks show the greatest rigidity [34J. 
Microtubule networks can withstand very high strain, and the greatest deformability. 
This is consistent with their role as structural support of the actin filaments |34j . 
Vimentin networks, which mostly make up the intermediate filaments, tend to be 
less rigid at low shear strain, but they harden at high strains. These responses make 
them ideal for the support of nucleus [31] . Because of this, we expect very significant 
differences in the responses of the actin cytoskeleton and the nucleus of the cell. 
Based on the above, we choose a simple network to model an endothelial cell. This 
model is shown in Figure 14.51 The flow sensor, body 1, is attached to the actin 
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cytoskeleton which is represented by bodies 2 and 3 in parallel. The actin cytoskeleton 
then attaches to the nucleus, body 4 in the diagram. 



1. 




Figure 4.5: Model I. of an endothelial cell. 



Next, we can elaborate on our initial diagram and add the connections between the 
nucleus and the attachments to the substrate. Part of this second model is the same as 
the first network, but now the nucleus (body 4) is further connected to actin bundles 
(bodies 5 and 6) that end at transmembrane proteins (body 7). 
1. 




Figure 4.6: Model II. of an endothelial cell. 



4.2.4 Parameter values 

Appropriate parameter values must be chosen for all of the bodies. The parameter 
values for the actin filaments are taken from Sato et al. [HI] who measure viscoelastic 
properties of endothelial cells with a micropipette technique. Guilak et al. give the 
parameter values for the nucleus based on a study also with a micropipette aspira- 
tion, and they conclude that the nucleus is about 3-4 stiffer and approximately twice 
as viscous as the cytoplasm [21]. Finally, the parameter values for transmembrane 
proteins is found by Bausch et al. ;6] who use the novel technique of magnetic bead 
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microrheometry. All the parameter values are summarized in Table 14.11 





m ( Pa s) 


Moi (Pa) 


Mn (Pa) 


Ref. 


Actin filaments 


5000 


50 


100 




Nucleus 


10 000 


200 


400 




Transmembrane proteins 


7.5 


100 


200 





Table 4.1: Parameter values for the endothelial cell models 



The parameter values for transmembrane proteins were given in the original paper in 
units of Pa m for the spring constants and Pa s m for the dashpot viscosity, and had 
to be converted to Pa and Pa s, respectively. These calculations are given below. 
Sato et al. give the formula for the deformation as 

L{t) = M1 e-*). 

TTMOI Mul + Mil 

The dimensions are as follows: [a] = m, [Ap] = Pa and [Moi>Mn] — P & - The formula 
for the deformation in Bausch et al. is 



L(t) = — (1 



Mil 



-e t ) 



Moi Moi + Mil 

with dimensions [F] = N, [moijMii] = Pa m = N/m. In order to compare the spring 
constants given in the two papers, we must have them in the same dimensions. First 
we will find what the applied force is in Sato et al., then we use the expression for 
the initial deformation and final deformation to determine the spring constants. The 
force is given by F = Apna 2 ~ 2500 pN. The initial deformation, 



L (t) = —(1 - 



Mn 



Moi Moi + Mn Moi + Mn 
and the deformation at steady state is 

F 

L s (t) = — . 
Moi 

Using the data in Sato et al. this gives us /z i = 6.35 x 10~ 4 Pa m and fin = 9.38 x 10~ 4 
Pa m. The same calculations with the data in Bausch et al. leads to /zqi = 1-25 x 10~ 3 
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Pa m and /jl u = 1.61 x 10 3 Pa m. Finally, we must obtain the viscosity of the dashpot, 

m = — : — • 

A*oi + A*n 

In Sato, rji = 4.125 x 10~ 2 Pa m s and in Bausch rji = 6.33 x 1CT 5 Pa m s. Now we can 
compare the parameter values, /xoi and //02 ior the nucleus is approximately twice the 
value of /i i and /i n , respectively, for actin filaments, and the viscosity in the nucleus 
is approximately 1.5 x 1CT 3 the dashpot viscosity of actin filaments. Based on this 
we arrive at: 

1*0! = 2(50) = 100 Pa 
/in = 2(100) = 200 Pa 
T] X = (1.5 x 10~ 3 )(5000) = 7.5 Pas. 
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4.3 Results 

This section contains two sets of numerical simulations. The first set of simulations 
examines the relationships between the parameters of the Kelvin bodies and the 
deformation, u(t) and the force splitting, a(t)F(t) in a the two-body problem. Both 
steady and oscillatory flow are considered. 

The second set of simulations takes the four-body and the seven-body model networks 
with realistic parameter values and finds the temporal dynamics of the deformation, 
and the force splitting coefficients. 

The differential equations are solved with a four-stage fourth-order Runge-Kutta 
method. The Matlab code used to solve the equations is presented in Appendix 
O The time step chosen for the simulations is h = 0.1. This method is of order 4, so 
the error is 0(h 5 ) = 0.00001. 

4.3.1 Parameter sensitivity analysis 

We investigate the behavior of two Kelvin bodies coupled in parallel. We focus on 
three questions: 

• temporal dynamics of the system with one parameter value changed over three 
orders of magnitude; 

• the steady state behavior of the system as a parameter value is changed over 
several orders of magnitude; 

• the behavior of the system when all three parameters are changed over three 
orders of magnitude 

This analysis allows us to understand how the parameter values determine the vis- 
coelastic properties of a material which gives us intuition into the behavior of the 
more complicated model networks involving materials of different properties. 
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In all of the figures depicting the deformation u(t) and the force splitting, a(t)F(t), 
the parameters for body one are set to baseline levels, /xqi = 50, /in = 100 and 
7711 = 5000. With the exception of Figures f4.3QII4.3cll in which all three parameters of 
body 2 are changed, in all other figures of this section two parameters of body two are 
kept constant, and only one parameter is changed. When the steady state behavior 
of the two Kelvin bodies is investigated, the following values are used: /io2 — £*i2 — 
1, 2.5, 5, 7.5, 10, 25, 50, 75, 100, 250, 500, 750 and 1000 and r] 12 = 1, 2, 5, 10, 20, 50, 
100, 250, 500, 1000, 2500, 5000, 10000, 25000, 50000 and 100 000. The baseline value 
for r)i2 is larger than the other two constants, and this is why we consider much larger 
values when examining the steady state behavior. The figures show which parameter 
is altered, and the values used in the simulations. The time steps of the simulations 
are h = 0.1, so one second of time is depicted by ten time steps on the graphs. 




0.6 0.8 1 1.2 

timesteps, h, h=0.1 



x 10 



Figure 4.7: Dependence of deformation on /xo2- Steady flow. 
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time, t 

Figure 4.8: Oscillations in the deformation for a particular value of /io2- The circles show the only 
points which are displayed in the subsequent figures displaying the dependence of deformation on 
parameter values. Oscillatory flow. 

Figures l4.7H4.10l show the deformation of the two body model as one of the spring 
constants, /i 2 changes. (/i 2 is the constant that characterizes the isolated spring.) 
Figure l4~7l shows the deformation in steady flow for three values of the spring constant 
/io2 =5, 50, 500. The spring becomes stiffer as /io2 increases, and this results in 
reducing the deformation, because the spring is more difficult to stretch. The steady 
state is reached very quickly for /i 2 large. 

Figure 14.81 illustrates the treatment of oscillatory flow. Because the frequency of 
oscillations is high, we do not want to display all the oscillations, only the peak values 
that the deformation reaches. These values are marked with an V on this figure, and 
in subsequent graphs of oscillatory flow only these are displayed. We also note in this 
figure that the steady state value in oscillatory flow is obtained within the first few 
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x 10"' 



1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 

time, t 



Figure 4.9: Dependence of deformation on /j, 2- Oscillatory. 



Steady state values of u(t) vs. n 




Figure 4.10: Dependence of steady state deformation on /Uq2- 
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oscillations (2-3) seconds. This observation is confirmed by figure I4~U1 which depicts 
the time evolution of the deformation for oscillatory flow as the spring constant /^ 2 
is varied. The steady state of the deformation is reached almost immediately, and it 
is a very small value (j 0.004), even if the spring constant is small /io2 = 5. Stiffer 
springs result in even smaller deformations, as expected. 

We can compare how the steady state of deformation changes with the spring constant 
/io2 in steady flow, and the steady state of the peak deformation in oscillatory flow in 
Figure I4~. 101 Each peak value is obtained after 2 x 10 4 time steps, which corresponds to 
2000 time units (seconds). Clearly, much larger deformation is produced in steady flow 
than in oscillatory flow, then, as the spring becomes increasingly stiff, the deformation 
tends to zero regardless of the flow. 

u vs t for various values of n 12 , F=F 

0.01 

0.009 

0.008 

0.007 

c 0.006 
o 

E 

£ 0.005 
o 

^0.004 
0.003 
0.002 
0.001 



0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 

timesteps, h, h=0.1 x1 g* 

Figure 4.11: Dependence of deformation on /ii2. Steady flow. 

The next three figures, Figures I4.11M.131 show how the deformation of a cell changes 
as the other spring constant of body two, /112 changes. (This spring constant charac- 
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. x 1 0"' 



u vs t for different values of n 12 , F=F cos(a) t) 



: 3 - 



H 12 =10 



H 12 =1000 



1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 

time, t 



Figure 4.12: Dependence of deformation on ^i 2 . Oscillatory flow. 



0.01 



Steady state values of u(t) vs. u. r 

t-¥ x ■ ■ >' ■ ' ¥ 




Figure 4.13: Dependence of steady state deformation on ^i 2 . 
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terizes the spring next to the dashpot in the second Kelvin body.) Figure 13. 1 II shows 
that although the initial deformation is different for the three values of /ii 2 (as one 
expects, because the initial condition depends on ^12), the steady state obtained is 
the same for the three values. Similarly, in Figure \%.12\ we see that in oscillatory flow 
the steady state is obtained virtually immediately again (in 2-3 seconds), and that 
the deformation for all three values of /ii 2 is very small, less than 0.005. 
In Figure 14.131 most steady state values are obtained after 2000 seconds again. In 
oscillatory flow, as the spring becomes stiffer, the deformation decreases, just like it 
did in Figure 14.101 However, in steady flow, the deformation stays almost the same 
regardless of the spring constant. This is due to more of the force concentrating on 
the dashpot if the spring is very pliable. In this case the steady state value is obtained 
slower too, in (3000 seconds) which also suggests that the force is concentrated on the 
dashpot which creeps to the steady deformation slower. Just like before, the overall 
deformation is always smaller in oscillatory flow, because the oscillatory flow only 
produces a large force periodically and not continuously. Also, there is a force of 
equal magnitude but opposite direction acting on the dashpot periodically, therefore 
the dashpot is unable to extend fully. (This claim is verified later in simulations where 
the frequency of oscillations is changed, allowing the dashpot more time to deform.) 
When the spring next to the dashpot is very pliable, i.e for small values of /ii2, much 
of the force is allowed to act on the dashpot, therefore the overall deformation is 
larger when the /ii 2 is changed ( Figure l4~T3*|) than if /z 2 is changed (Figure l4.10|) . 
Now we turn to examining how the dashpot viscosity influences the behavior of 
the system. Figures 14.141 and 14.151 depict the time evolution of the deformation for 
?7i2=500, 5000 and 50000. In steady flow the steady state of the deformation is the 
same regardless of the value of 7/12, but for large dashpot viscosities the steady state 
takes longer to obtain, approximately 5000 seconds. Figure 14.151 shows the steady 
state of the deformation in oscillatory flow. For the three values of 7712 shown here, 



CHAPTER 4. 



MATHEMATICAL MODELS IN BIOLOGY 



148 



u vs t for various values of ii 12 , F=F 




0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 

timesteps, h, h=0.1 10 i 



Figure 4.14: Dependence of deformation on 7712. Steady flow. 

the difference in deformations is negligible, on the order of 10~ 6 . The largest defor- 
mation is obtained when ry 12 is the smallest, and the bodies in this case deform very 
quickly, within the first 100 seconds. 

Figure 13.161 is created by running the simulations longer than in the previous figures, 
to time t — 1.2 x 10 4 seconds to ensure that the deformation reaches its steady state. In 
oscillatory flow, the steady state of the deformation is slightly larger for small dashpot 
viscosities, but after about 7712 = 100, the deformation is independent of the viscosity 
The interpretation of this is that when the viscosity is sufficiently low, the force is 
able to deform the dashpot quickly, then, as the forcing oscillates the deformation 
remains the same. For larger viscosities the initial deformation of the dashpot is 
negligible, and the later oscillations are again unable to change the deformation of 
the dashpot. As noted before, the deformation in steady flow is independent of the 
dashpot viscosity. 
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Figure 4.15: Dependence of deformation on 7712. Oscillatory flow. 
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Figure 4.16: Dependence of steady state deformation on 7712- 
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a vs t for various values of n Q2 , F=F Q 




Figure 4.17: Dependence of force splitting coefficient on /xq2- Steady flow. 



Now we can turn to investigating how the force splitting is effected by the parameter 
values. All the figures of the force splitting show the force acting on body one, and the 
force on the second body can be found from this by recalling that the sum of the two 
forces is one. We use the same values for /X02 as previously. The simulations in steady 
flow are shown in Figure [4". 171 As expected, if the spring is very pliable in body two, 
then most of the force will have to focus on body one. For two identical bodies (i.e. if 
fiQ2 = 50), the force splits equally between the two bodies. A very stiff spring in body 
two means that all the force has to concentrate here. In oscillatory flow we are only 
interested in the peak force acting on body one. Figure 14.181 shows the oscillations of 
the overall force, the force acting on body one, and the force splitting coefficient, a 
displayed on the same graph. Only the peak values of aF (which are identical to the 
value of a) are displayed subsequently. As before, the peak values of the force do not 
change with time, as shown in Figure 14". 191 In the oscillatory flow, just like in steady 
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aF,F and a vs t for F=F cos (to t) for identical bodies 
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Figure 4.18: Oscillations of a(t)F(t), F(t) for identical bodies. a(t) also displayed. V marks the 
points which are used to display peak oscillations of a(t) in the following figures. Oscillatory flow. 

flow, a pliable spring in body two leads to more of the force concentrating on body 
one, in identical bodies the forces split evenly, and stiff springs require more force to 
deform. 

Figure I4.2UI depicts the steady state values of force splitting for a range of values of 
/i 2- The curves for steady flow and oscillatory flow intersect when fi 02 =50, because 
this is the baseline value at which the two bodies are identical, therefore the forces 
split evenly regardless of the flow. It is interesting to note that in steady flow the 
force splitting is more extreme than in oscillatory flow. More specifically, a small /io2 
(pliable spring in body two) in steady flow allows much more of the force to act on 
body one than in oscillatory flow, but a large value of //02 i n steady flow leads to a 
smaller force on body one than in oscillatory flow. The interpretation of this is that 
in oscillatory flow the dashpot offers a constant resistance, thus much of the overall 
force is always trying to stretch the dashpot. 
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Figure 4.19: Dependence of a(t) on fi 2- Graph shows peak values of a(t)F(t). Oscillatory flow. 



Steady state values of aF(t) vs n 




Figure 4.20: Dependence of steady state force splitting on /i 02 . 
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Figure 4.21: Dependence of force splitting coefficient on ^12. Steady flow. 

Next, we examine the dependence of the force splitting on the other spring constant, 
fiu- Figure l4*.21l shows the force splitting in steady flow. The initial condition of the 
force in the bodies depends on /ii2, but clearly, the steady state of the force tends to 
the same value, 0.5. Regardless of the stiffness of the spring in body two, eventually 
the forces acting on the two bodies become the same. When the force acts on a pliable 
spring, more force goes to the dashpot, and when the spring is stiff, more of the force 
goes to it. The mediating effects of the dashpot lead the equal force splitting between 
the two bodies. 

In oscillatory flow, shown in Figure 14.221 the mediating effects of the dashpot are 
smaller, so there is a larger difference between force splitting between bodies one and 
two for different values of the spring constant, fi\2- As before, the peak force acting 
on body one reaches its steady state very quickly, so we see constants. For small 
values of /112 the more of the force acts on body one, and for stiff springs more force 
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aF vs t for different values of \i , F=F cos(co t) 
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Figure 4.22: Dependence of a{t) on /ii 2 . Graph on shows peak values of a(t)F(t). Oscillatory flow. 

Steady state values of aF(t) vs n 12 




Figure 4.23: Dependence of steady state force splitting on /x 12 . 
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is concentrated on body two. The comparison of steady states in oscillatory flow and 
steady flow, as depicted by Figure 14.231 brings no significant new information. 



a vs t for various values of ii 12 , F=F 



0.7 - 



£ 0.6 



S 0.5 



Tl =50000 



0.5 



1.5 2 2.5 3 

timesteps, h, h=0.1 



3.5 



4.5 



X 10 



Figure 4.24: Dependence of force splitting coefficient on 7/12- Steady flow. 



Now we turn to looking at the dependence of force splitting on the dashpot viscosity. 
Longer simulation times are necessary again to obtain the steady state values. In 
steady flow, shown in Figure 14.241 small dashpot viscosity in body two leads to a 
transient increase of the force acting on body one, but at the steady state the force 
splits equally between the bodies. The transient increase is due to more force needing 
to deform the dashpot of body one. When the dashpot viscosity of body two is large, 
there is a transient decrease in the force acting on body one. The reason for this is 
similar to our argument above: now more force is necessary to deform the dashpot 
of body two. It is also clear from the figure that for large dashpot viscosities the 
transient time is longer. This is consistent with our previous observations on the role 
of the dashpot viscosity. 
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Figure 4.25: Dependence of a(t) on rji 2 . Graph shows peak values of a(t)F(t). Oscillatory flow. 

Steady state values of a(t)F(t) vs. r\ i2 
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Figure 4.26: Dependence of steady state force splitting on fy 12 . 
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Dependence of peak deformation on frequency 
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Figure 4.27: Dependence of peak steady state deformation on the frequency. Peak steady state 
deformation of oscillatory flow is divided by the steady state value of steady flow. 

In oscillatory flow, Figure 14.251 the force split is almost exactly equal independently 
of the dashpot viscosity Small viscosities lead to a quicker and slightly larger de- 
formation. Figure I4.26I shows how the steady state values of the force in body one 
change with the dashpot viscosity in steady and in oscillatory flow. The most notable 
observation is that in oscillatory flow the forces do not split equally. This is the result 
of the initial quick dashpot relaxation of body two (when rji 2 is small) that leads 
to more force being necessary for body one. Similarly, if the dashpot viscosity of 
the second body is large, then the initial dashpot relaxation occurs in the first body, 
therefore more force will always be applied to the second body. 
So far all of our numerical simulations in oscillatory flow used the frequency / = 
1Hz = 2tt rad/sec. This is a physiological value which corresponds to the frequency 
at which the heart pumps the blood through the blood vessels, however, this frequency 
may change during exercise, or due to pathologies. This raises the question of how 
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Dependence of peak deformation on frequency and r\ 
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Figure 4.28: Dependence of peak steady state deformation on the frequency and 7712. Peak steady 
state deformation of oscillatory flow is divided by the steady state value of steady flow. 

the frequency of oscillations may change our model, in particular how the frequency 
effects the deformation and force splitting of the Kelvin bodies. 
Figures f4 . 2 71 - I4~2"§1 show the results of simulations where the frequency of oscillations 
is altered. In these plots we used the values / = 10~ 4 , 2 x 10~ 3 , 3 x 10 -3 , 4 x 10~ 3 , 5 x 
1(T 3 , 1(T 2 , 2 x 1(T 2 , 5 x lO- 2 , 7 x 1(T 2 , 10 _1 , 2.5 x 10 _1 , 7.5 x 10 _1 , 1, 10 Hz. All other 
parameter values for body 1 and body 2 are baseline values. The peak deformation 
(and the peak force splitting) for each of the simulations is divided by the deformation 
(and the force splitting) for steady flow. 

Figure 14.271 shows that if the frequency of oscillations is very small, then the defor- 
mation in oscillatory flow is essentially the same as in steady flow. As the frequency 
increases, the deformation decreases, because the dashpot is unable to fully deform 
once the oscillations become fast enough. This is the result of the sign of the force 
changing quickly, therefore the force not acting for sufficiently long periods of time 
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Figure 4.29: Dependence of peak steady state force splitting on the frequency. Peak steady state 
force splitting of oscillatory flow is divided by the steady state value of steady flow. 

to fully stretch the dashpot. Once the frequency reaches a critical value, f crit ~ 10~ 2 , 
the deformation is independent of the frequency of oscillations. This signals the fre- 
quency at which the deformation is entirely due to the springs. (Deformation of the 
springs is instantaneous when force is applied.) 

Just like in the one Kelvin body case, f cr u = 10~ 2 two orders of magnitude 
below the physiological value. This implies that endothelial cells exposed to purely 
oscillatory flow will only deform to a small fraction (approximately one third) of 
deformation possible in steady flow. If there is a threshold value of deformation that 
permits cells to align and go through other significant physiological changes, it is very 
likely that purely oscillatory flow will not be able to reproduce these effects. Also, it 
is interesting to note that the deformation for frequencies / > f cr u is the same for 
one Kelvin body as for two. This implies that any number of Kelvin bodies coupled 
in parallel will not deform more than this value, therefore a model consisting of n 



CHAPTER 4. MATHEMATICAL MODELS IN BIOLOGY 160 

u vs t, each parameter changed by a factor of 1 0, F=F 

0.02 
0.018 
0.016 
0.014 

c 

O 

| 0.012 

jo 
"o 

!? 0.01 
o 
c 

<D 

^ 0.008 
0.006 
0.004 
0.002 


0.5 1 1.5 2 2.5 3 

timesteps, h, h=0.1 x1 q* 

Figure 4.30: Deformation when every parameter of body 2 is changed by a factor of 10. Steady flow. 

Kelvin bodies in parallel would retain the feature that f crit is well below physiological 
values. 

We have heuristically argued before that for all values of viscosity 7712 > 100 the 
dashpot is not able to react in oscillatory flow, because the direction of the force 
changes very quickly, and there is not enough time for the dashpot to deform. This 
suggests that the relationship between frequency and deformation has to be exam- 
ined. Figure 14". 281 shows these simulations. We observed in the previous plot, (Figure 
14.27(1 that as the frequency increases, the overall deformation decreases to a constant 
value. This can be explained, if very small frequency oscillations allowed sufficient 
amount of time for the dashpot to react, but once the frequency increased to about 
/ = 0.01 Hz, the dashpot did not have time to respond at all. Beyond this frequency 
all the deformation is due to the springs. Our figure confirms that by making the 
dashpot very inviscid (decreasing r/ 12 below 100) allows the dashpot to react much 
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Figure 4.31: Deformation when every parameter of body 2 is changed by a factor of 10. Oscillatory 
flow. 

more quickly, and we see larger deformations for given frequencies. When the fre- 
quency is very low, then the deformation is similar to the deformation in steady flow 
for any value of the viscosity, therefore the dimensionless value (peak deformation 
in oscillatory flow divided by the deformation in steady flow) is near 1. When the 
frequency is very large, then the deformation is always a constant around 0.333 in- 
dependently of the frequency or the viscosity. Between these two regimes there is a 
range of frequencies for which increasing the frequency means decreasing the defor- 
mation. Interestingly, for small values of the dashpot viscosity the normalized peak 
deformation is bi-sigmoidal. 

Figure 14.291 describes the dependence of the peak force splitting on the frequency 
of oscillations. As before, the largest amplitude at the steady state is taken for the 
appropriate value of the frequency, and this amplitude is divided by the steady state 
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aF vs t, each parameter changed by a factor of 10, F=F 
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Figure 4.32: Force splitting when every parameter of body 2 is changed by a factor of 10. Steady 
flow. 

value in oscillatory flow. As we can see, the frequency of oscillations does not change 
the fact that the peak force split in oscillatory flow is always the same as the steady 
state force split in steady flow. 

We must also examine the case that all parameters in the second body are changed, 
because contradictory predictions could be made based on changing individual pa- 
rameters only. Figures r4.30N4. 3^1 reveal that regardless of the type of flow, the defor- 
mation decreases if the parameters are increased, and the deformation is increased if 
all the parameters are decreased. Similarly, in either type of flow the force in body 
one decreases if the parameters are increased in body two, and the force increases in 
body 1, if the parameters in the other body are decreased. 
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a vs t, parameters changed by a factor of 10, F=F cos(co t) 
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Figure 4.33: Force splitting when every parameter of body 2 is changed by a factor of 10. Oscillatory 
flow. 

4.3.2 Network simulations 

We represent endothelial cells as a network of viscoelastic bodies. Each part of the cell 
we model, namely, the transmembrane proteins, flow sensors, cytoskeletal elements 
and the nucleus are thought of as Kelvin bodies with different parameter values for 
the spring constant and the dashpot viscosity. The baseline values which we use for 
actin filaments are: /J02 = 50 Pa, fi\2 = 100 Pa and 7/12 = 5000 Pa s. For the nucleus 
the spring constants are four times the baseline values and the viscosity of twice the 
baseline: /io2 = 200 Pa, fiu = 400 Pa and 7712 = 10000 Pa s. These estimates are 
based on experimental measurements by Guilak et al. [21] . The parameter values for 
the transmembrane proteins (flow sensors, ion channels, attachments to the substrate) 
based on |6| are: /i 2 = 100 Pa, fii 2 = 200 Pa and r] 12 = 7.5 Pa s. The parameter 
values are calculated in Section T4.2.4I and summarized in Table I4~T1 
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Figure 4.34: Deformation of different parts of a simple endothelial cell in steady and oscillatory. 



Figure I4.34I shows numerical simulations for the simple four-body model of the en- 
dothelial cell fFigure l4~o"|) in steady and in oscillatory flow. In steady flow transmem- 
brane proteins deform immediately, followed by the deformation of the nucleus, then 
the deformation of the actin filaments. This behavior is consistent with transmem- 
brane proteins having the smallest, and actin filaments having the largest dashpot 
viscosities. Transmembrane proteins and actin filaments reach the overall deforma- 
tion whereas the nucleus only deforms slightly. Because our model is linear, the 
overall cell deformation is the sum of the deformation of the components. The time 
constants are consistent with experimental data in which transmembrane proteins 
such as flow sensors and ion channels respond to flow very quickly, on the order of 
seconds, and cytoskeletal reorganization is the slowest, in fact, it happens on the time 
scale of many hours. 

Next, we examine the same four-body model in oscillatory flow. Only peak values 
of the deformation are plotted, as before. The overall deformation is much smaller 
here than in steady flow, and the steady state of deformation is attained immediately 
(within 2-3 seconds). The nucleus deforms the least amount again, and here the 
largest deformation is that of the transmembrane protein. 
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Because in our four-body model the nucleus and the flow sensor are both modeled 
as a single Kelvin body, the force acting on each of them is a constant, F = 1. The 
same parameter values characterize the two actin filaments which are modeled as two 
Kelvin bodies connected in parallel, therefore the forces acting on them are also equal, 
aF = (1 — aF) = 0.5 This is true for steady as well as in oscillatory flow. 
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Figure 4.35: Deformation of different parts of a simple endothelial cell in steady and in oscillatory 
flow. 

The final graph, Figure I4.H5I and shows the response of the seven Kelvin-body model 
(depicted in Figure E~T)|) in steady and oscillatory flow. The results are very similar to 
the simulations of the four-body model. The forces for the bodies are all constants: 
for the transmembrane proteins and nucleus (which are modeled as single Kelvin 
bodies) the force is equal to one, and for acting filaments, modeled as Kelvin bodies 
in parallel the force is split equally. 
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4.4 Discussion 

Our numerical results included a parameter sensitivity analysis and simulations for 
two simple model networks of endothelial cells. The results of the simulations have 
been interpreted so far only in terms of the elements of the Kelvin body, but it has not 
been discussed what their implications are for endothelial cells. This section describes 
the conclusions we can draw regarding endothelial cell behavior from the numerical 
results. We also briefly mention further work that can be done to improve the current 
model. 

The most notable difference between the effect of steady and oscillatory flow is that in 
oscillatory flow deformation tends to be much smaller than in steady flow for the same 
viscoelastic materials, as it is depicted in Figures ^. 101 14.151 and I^TBl The importance 
of the large deformation difference in the two types of flow is clear when we recall 
experimental results in which endothelial cells exhibit elongation and certain other 
biological responses (such as activation of flow- sensitive ion channels and mobilization 
of intracellular calcium) depending on the specific form of the shear stress cells are 
exposed to. Dewey et al. original source: [T7]) demonstrate that endothelial 

cells exhibit some morphological responses when exposed to large shear stress, for 
example in steady or pulsatile flow, but not when their exposure to shear stress is 
below a certain level, for example, in oscillatory flow. 

According to our model, the forces acting on each part of an endothelial cell, regardless 
of its viscoelastic properties elicits a much smaller response in oscillatory flow. This 
can be particularly important for a flow sensor, as discussed by Barakat which 
in oscillatory flow may not be stimulated sufficiently to initiate a signaling cascade 
for downstream responses to the flow. Our model shows, that even given an signal 
from the flow sensor, the deformation of materials whose viscoelastic properties are 
consistent with those of the nucleus and cytoskeletal elements would be much smaller 
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in oscillatory flow. 

Another observation we can make (still based Figures 14. 1(J| 14.131 and 14.16)) is that 
the viscoelastic materials which display the largest difference in deformation can be 
characterized by small spring constants but a viscosity coefficient which is at least 
7] « 10 2 . This characterization would allow us to estimate parameter values for 
microtubules and intermediate filaments for future simulations, because we have ex- 
perimental observations of the relative behavior of actin filaments, microtubules and 
vimentin 34 . 

We have investigated how the frequency of oscillations effects the deformation. Our 
simulations show that within physiological conditions (even accounting for conditions 
where the frequency of the flow changes one order of magnitude) the normalized de- 
formation is small, unless the viscosity coefficient becomes very small. This implies 
that materials whose viscosity is small, for example various transmembrane proteins, 
such as attachments to the substratum or cell-cell adhesions, are less able to differ- 
entiate between oscillatory and steady flow than the cytoskeleton or the nucleus. (As 
discussed above, the cytoskeleton and the nucleus respond very distinctively to steady 
and to oscillatory flow.) 

Now let us examine the network simulations of endothelial cells. Both networks 
hugely oversimplify the complex connections between endothelial cells, therefore no 
detailed conclusions should be drawn from our results, rather we ought to make 
simple qualitative observations. More sophisticated networks need to be created in 
the future, and the main significance of our current simulations is to demonstrate the 
sort of models we are able to create with our elements now (namely n Kelvin bodies 
in parallel, coupled in series with single Kelvin bodies). 

Qualitatively, the time scale of deformations as predicted by our model fits with ex- 
perimental observations. Instantaneous response can be seen from transmembrane 
proteins (flow sensors, focal adhesions, cell-cell adhesions), followed by the deforma- 
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tion of the nucleus, and finally, the changes in the cytoskeleton. Flow sensor and ion 
channel deformation can occur within seconds whereas changes in the gene expression 
takes hours, and the cytoskeletal re-organization takes place over the span of about a 
day. These time scales do not compare well with the prediction of our model, which 
does predict flow sensor and ion channel deformations within seconds, but it also 
predicts the response of the nucleus and the actin filaments to be on the order of 
minutes. Clearly, our model must be modified, and in order to obtain quantitatively 
accurate results, the effects of the biochemical signaling pathways must be considered 
as well. 

Our model is the first step in creating more realistic model networks of endothelial 
cells. The current project can be extended in a number of ways to lead to better 
approximations of the cellular response to flow. The extension promising the most 
extensive changes in the results is the assumption that the parameters describing 
the Kelvin bodies can alter in time too. This new assumption would turn our linear 
system of equations to a nonlinear system with possibly much more complicated 
dynamics. Other changes to incorporate in the model are: formulating the equations 
for new networks based on combinations of bodies which our system does not describe, 
assuming that the forcing function is applied gradually (there is a "ramp" for the 
force). Allowing the force to act gradually on the bodies could be compared to 
experimental results in which cell responses are altered by applying shearing forces 
instantaneously or gradually. Although even the understanding of viscoelastic Kelvin 
bodies in parallel lead to new insights on the hypothesized mechano-transduction of 
flow induced shear stress in endothelial cells, further improvements and extensions 
are necessary to clarify the role of this mechanism. 
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Appendix A 

Receptor model 

An important question is whether there is biological evidence of the signal transduc- 
tion pathway allowing the existence of the proposed turning rates that depend on the 
oxygen concentration and the spatial gradient of oxygen. In [22], Taylor and Johnson 
propose a model that aims to explain experiments in which receptors are rewired to 
produce inverse responses. Using their model for a receptor, it is easy to explain how 
the turning rates of our model could be a result of a simple two-state receptor. 




Figure A.l: The four-helix bundle of Tar receptor. (Figure from Taylor and Zhulin, [73].) 
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It is known that part of the Tar (and Tsr) receptors is a four-helix bundle outside the 
cell, and as ligands bind to the receptor, the bundle goes through a conformational 
change. By representing the receptor as two independently moving parts, we have a 
simplified model ('piston model') of the receptor. Each part of the receptor is able to 
move in the z direction as a function of the proton motive force, PMF. One part of 
the receptor is able to change its position fast in response to PMF, while the other 
part reacts slower. When the two parts of the receptor in the model are at the same 
position, they lock (representing the conformation change in the real receptor), and 
this initiates tumbling. When the two parts of the receptor are not in the same 
position, the cell swims smoothly. 




Figure A. 2: Piston model of a receptor. (Figure based on Taylor and Johnson, |72|.1 
In ground state, the two parts are displaced from one another. Thermal energy is 
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able to randomly move the parts, which causes tumbling in unstimulated cells if the 
parts are moved together, and it promotes smooth swimming if the parts are moved 
further apart. Attractant binding to the receptor also moves the fast and slow parts 
further, which makes smooth swimming more likely. Adaptation returns the receptor 
to the ground state. 

In our mathematical model such a receptor could be responsible for the same oxygen 
concentrations triggering different signals. When a bacterium is outside the optimal 
oxygen concentration, the cell is in ground state; therefore, it has its baseline turning 
frequency. As the cell enters the favorable oxygen concentration, the fast-changing 
part moves; thus, the two parts get further from each other, and the cell swims 
smoothly inside the band. By the time the cell gets across the band, the slow-changing 
part moves too, and the two parts lock which causes tumbling. The tumbling turns the 
cell back into the favorable concentration, and the fast-moving part changes position 
again, leading to smooth swimming. 

Consider zf. z coordinate of fast-changing part, z s : z coordinate of slow-changing 
part. 



In these equations i} and z s are the equilibrium positions of the fast-moving and 
slow-moving parts, respectively. In ground state, the equilibrium positions are not 
the same, as discussed above, zf and zf^ are constants, p is the proton motive force 
(which is proportional to the oxygen concentration). 



We can assume that the PMF is a linear function of time, since the fully developed 
oxygen gradient is linear in space, and it is observed by the bacteria swimming through 
it as oxygen changing linearly in time. 



Zf — Zf + Cip 
;(°) + Cl p = zf - Az + Cip 




Zs = Z. 



(A.2) 



p = ±kt + c 
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fast (signaling domain/ 
critical residue) 



slow (catalytic domain/ 
active site) 




PMF 




PMF 



Figure A. 3: Fast and slow states of a receptor. The top figure shows how the receptor switches from 
high turning frequency to low turning frequency and vice versa. The bottom figure shows a different 
representation of the same process. (Bottom portion of the figure based on Taylor and Johnson, 

We can assume that the fast-moving part is just a function of the proton motive 
force (i.e. its position changes immediately as it detects changes in the proton motive 
force.) We assume that the slow-moving part changes its position after a delay. This 
results in the following two equations: 



z f « Zf(j)(t)) 

dz s 1 

~7T = -(Zs{p(t)) ~Z a ) 

at t 



(A.3) 
(A.4) 



The first equation is an algebraic equation, and we can solve it by substituting the 
expression for the proton motive force into (jA.3|) . This gives the equation: 



Zf = [z^p + CqCi] ± C\kt 



This shows that the coordinate of the fast part will be increasing as it moves up the 
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concentration gradient, and it will be decreasing as the cell swims down the gradient. 
Now we can solve (|A.4|) . 

dz s -z s 4 - Az + cip 
dt r r 
We have to solve the homogeneous equation on the left-hand side, then find the 
particular solution, z sp = A + Bt, using the right-hand side of the equation: 

dz s -Zs _ _ AZ + C l ( ±H + C °) 

dt T T 

Solving the homogeneous equation we get: z sh = z e~^ r - In the particular solution, 
we get A = z^ + c\Cq =f c\kr — Az; B = ±cik. This gives the full particular solution, 
z sp = z^p + c\Cq =f C\kr — Az ± C\kt. By making the appropriate substitution, this 
results in the full solution, 

z s = Zf + [TCikr - Az] + z e~ 

If the adaptation time, r, is small compared to the run (which is our assumption in 
the model), then the last term is negligible. Now we can see how the receptor will 
work by examining the position of the fast- and slow-moving parts. The position of 
the slow-moving part, z s is behind the position of the fast part. To see the dynamics 
of the two parts, the solution to z s can also be re-written as: 

z s = z^ + CiCq — Az ± c\k{t — r) + z e~ 

As long as the time, t is less than the adaptation time r, the slow part will be 
behind the fast one, but after the time reaches r, the sign of C\k(t — r) changes, and 
the difference in the position of the fast and slow part decreases. This difference is 
maintained until the cell turns down the gradient, at which point the z-coordinate of 
the fast moving part will change quickly again, followed by the z-coordinate of the 
slow moving part. 

This simple model of the receptor based on proton motive force is sufficient to explain 
the turning frequencies assumed in the full mathematical model of aerotaxis. 
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Appendix B 



Analytical calculations 



B.l Asymptotic approximation 

We want to solve the system of equations 

dM 



= m + \(-k a (h)M + k d A) 



dt 

with initial conditions: 



dt 
dA 

' -rA + X(k a (h)M-k d A) 



rk a (l ) r 

On the fast time scale, i.e. when r = At, 

dM dA 
dr dr 

which implies A + M = C for some constant, C. We substitute this expression into 
the original system of equations to get: 

dM 



(h -k a (h)M + k d (C-M) 
dA 

— = k a (h)(C-A)-k d A 
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The initial conditions on the fast time scale are the same as the initial conditions of 
the original equations, 

mk d 



M = M(r = 0) 



rk a (l ) 



A s = A(r = 0) = - 
r 



The solution to these equations is given by: 



M(r) = (M -^-)e-^ + ^ 

Tf Tf 

Where we defined the fast time scale to be 

rf = k a (li) + k d . 

Now we use the fact that the sum of A and M is always a constant, in particular, 
A(0) + M(0) = C. 

m ^ k d 



Two new constants are useful: 

k a (h)C ml + (k d /k a (l )) 



A 1 
Mi = 



r f r 1 + (k d /k a (h)) 
k d C m k d 1 + (k d /k a (l )) 



r f r k a (h) 1 + (k d /k a (h)) 
Substituting the expressions for A\ and M\ into eqn. IB. II and IB. II we have obtain A 
and M on the fast time scale: 

A(t) = (A s - A 1 )e' rsT + A x (B.l) 
M(r) = (M - Mx)e~ r f T + M x (B.2) 
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We turn to the slow time scale, and we note that now A = ^4^M. Substituting this 
expression into the original equations gives us: 

dM 



dt 

k a (h) dM k a {l 



-r- 



k d dt k d 

By adding equations IB. 31 and IB. 41 we get: 

dM _ mk d rk a (h) 

dt k a (h) + k d k d + k a (h) 

M(t) = ce^ + 



m (B.3) 
i^M (B.4) 



(k a (li) + k d )r s 
M(t) = ce ( - rat) + M 2 (B.5) 



We have defined 

k a (h) 



r, = r- 



kd + k a (h) 
and 

mk d 



M 



2 



rk a (h 

Similarly, we have an equation for A: 



A (t) = ce(-**) + HI (B.6) 

k d r 



The solutions on the fast and slow times scale must match, so the limit of eqn IB.l 
as r — ► oo must be equal to eqn IB.6I at zero. Similarly, \\m. T ^, 00 Mj = M s (t = 0). 
Matching the solutions allows us to find the expression for c. 

k a (h)c 



m 
r 



k d 

k d / . k d \ kd 



+ A S = A X (B.7) 



The full solution is the sum of the fast and slow terms with the common limit [A\ 



and Mi) subtracted. 

A{t) = A S + (A s - A^e^ + {A 1 - A s )e~ rst (B.9) 
M(t) = M 2 + (M - Mi)e~ r / At + (Mi - M 2 )e~ rst (B.10) 
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B.2 Steady state solution 

We start with the four equations describing the two compartment model of the signal 
transduction system again. 

^± = m + X(-k a (h)M 1 + Mi) + h(M 2 - Mi) (B.ll) 
at 

HA, 

— - = -rA x + A(A; (/ 1 )M 1 - Mi) + h(A 2 - A x ) (B.12) 
at 

—l = m + X(-k a (l 2 )M 2 + k d A 2 ) - h(M 2 - M x ) (B.13) 

^ = - r A 2 + X(k a {l 2 )M 2 - k d A 2 ) - k 2 (A 2 - A,) (B.14) 
at 

Afk(O) = M 2 (0) = A 1 (0) = A 2 (0) = - (B.15) 

r k{t ) r 

At the steady state the left hand side of these equations is zero. We start finding the 
expressions for the steady state by adding equations IB.11I and IB. 131 In the following 
calculations we write A±, A 2 , Mi and M 2 instead of Au, A 2s , M\ s and M 2s . 

2m - Xk a (h)M 1 - Xk a (l 2 )M 2 + Xk d (Ai + A 2 ) = (B.16) 

2?77 

k a (l 1 )M l + k a (l 2 )M 2 = — + k d {A x + A 2 ) (B.17) 

Next we add equations IB. 121 and IB. 141 

- r{A x + A 2 ) + X{k a {h)M 1 + k a {l 2 )M 2 ) - Xk d (A 1 + A 2 ) = (B.18) 
(Ax + A 2 ){Xk d + r) = X(k a (l 1 )M 1 + k a (l 2 )M 2 ) (B.19) 

We substitute the expression for k a {l\)Mi + k a {l 2 )M 2 from lB.lTI 

A 1 + A 2 = — + k d (A! + A 2 , 

r + Xk d \ X 

2m Xk d / \ * \ / \ 

A > + A > = 7TW d + 7TWj A + A < } (&20) 

Solving the equation for A\ + A 2 we obtain 

A, + A 2 = — (B.21) 

r 
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This allows us to express A 2 (A 1 ) = ^ — A\. 

Now we return to IB. 171 to find an expression for k a (l\)M\ + k a (l 2 )M 2 explicitly. 



2/72 2771 

fc«(Zi)Mi + k a (l 2 )M 2 = — + k d — 

A r 

2m (r + Xk d \ fc (/i),, 
M 2 = — — — - — - - t^TT^I 



k a (h) V Ar / k a (l 2 ) 
Now we add equations IB. Ill and IB. 121 and similarly, add IB. 131 and IB. 141 

m - rAi + h(M 2 - M x ) + k 2 (A 2 -A 1 )=0 
m-rA 2 - ki(M 2 - Mi) - k 2 (A 2 - A x ) = 



By subtracting IB . 25l from IB~24l we arrive at 

Ai-A 2 



(M 2 - Mi) 



(B.22) 
(B.23) 

(B.24) 
(B.25) 

(B.26) 



r + 2k 2 

We want to express A± as function Mi, and this will allow us to find A 2 (Mi). In 
order to do this, we add IB. 211 and IB. 261 



A x = - + —±—(M 2 -M x ) 
r r + 2k 2 

Now we use IB.231 to find both A\ and A 2 as a function of M 1; so we have A\{Mx t 
A 2 {M X ) and M 2 (Af 1 ). 



a 1 = ™+ ki r 2m(r+Afcrf) -Mi) Ml _ Ml 



A, 



r r + 2A; 2 V A; a (/ 2 )Ar fc tt (fe)' 
fci (2m(r + \k d ) k a (h) 



rn 



Mi-M 



(B.27) 
(B.28) 



r r + 2/c 2 V k a (l 2 )Xr k a (l 2 )' 
We return to the equation IB. Ill with its right hand side set to zero, and substitute 
Ai(Mi), A 2 (Mi) and M 2 (Mi) from the equations lB~27l E2HI and lB~23l respectively, 
and solve the equation for Mi. 



777 + A 

kg(h) + k a (l 2 
k a (h) 



/ (i \ a/t , i f m , fcl r 2"7(r + A/c d ) 
- + fc d (- + [ fca(/a)Ar 



K(l 2 )Xr 



k a (k 
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mXk d 2m(r + Xk d ) 
m H h ky 



Xk a (h) + h 



k a (h)Xr 
k a (k) + k a (l 2 



Xk d 



(; 



r + 2k 2 
\k d 



+ 1 



1) 



Mi 



k a (h) K r + 2k 2 

If we simplify this expression, and substitute it back into IB. 271 IB. 281 and IB. 231 we 

arrive at the steady state solution: 

T\k\k 



At 



M 1 



Mo 



m 
r 
m 

A 2 = — 

r 



Xr 
mri 
Xr 



Xr 2 k p + hk s (r 2 + Xk a 
—T\k\k 



1 + 

1 + Xr 2 k p + k 1 k s (r 2 + Xk d ) 
r 2 {Xk a {l 2 ) + 2k l ) + 2Xk d k 1 

- r 2 [Xk p + kik s ] + Xk d kik s 
-r 2 (Xk a (l 1 ) + 2k 1 ) + 2Xk d k 1 

- r 2 [Xk p + kik s ] + Xk d kik s 



(B.29) 
(B.30) 
(B.31) 
(B.32) 



where we have defined: 



ri = r + Xk d , r 2 = r + 2k 2 
k = k a (h) - k a (l 2 ), k s = k a (li) + k a (l 2 ) 

k p = Ki^kaik) 

We want to verify that for k a {li) = k a {l 2 ) = k a we obtain the original steady state. 
This is clear for A\ and A 2 by inspection, but we need to simplify M\ and M 2 . We 
show the calculations for M\. 



Mi 



mri 
k a Xr 



r 2 (Xk a + 2k 1 )+2Xk d k l 
r 2 [Xk a + 2k 1 ] + 2Xk d k 1 . 

M = m(j + Xkd 



lim Mi 

A— >oo 



Xrk a 
mk d 



rka 



Previously we have shown that if k± = then the two compartments respond to 
stimulus as if they were not connected. In the main text we also mention that no 
important qualitative changes occur when k 2 = 0. In this case our system of equations 



APPENDIX. 



MATHEMATICAL MODELS IN BIOLOGY 



180 



becomes: 



m 

A 1 = — 

r 



1 + 



T\k\k 



\rk„ + k-ik s r\ 



A 2 = 



m 



1 



r\k\k 1 
Xrkp + k\k s ri\ 



M 1 
M 2 



Xr 



r(\k a (l 2 ) + 2fci) + 2Afc d fci 



dkik s J 



r[A/c p + fcifcg] + A/c, 
r(\k a {li) + 2fci) + 2Afc rf fci 



r[AA; p + fcifc a ] + Xk d kik s 
It is clear from the above equations that A 1 and A 2 depend on the ligand difference, so 
the system will respond to ligand gradients. The main text shows that the assumption 
k 2 ^> 1, on the other hand, results in cells where A\ — A 2 , so the cell cannot maintain 
an internal ligand gradient. 

We also want to find the steady state of the system for A 3> 1, and show that the 
qualitative behavior remains the same as in the A 0(1) case. We rearrange Mi 
and A\ to show decreasing powers of A, and note that similar rearrangements can be 
made for M 2 and A 2 . 

m X 2 a + Xb + c 
1 r X 2 d + Se 
a = k a (l 2 )k d (r + 2k 2 ) + 2k x k\ 

d=(r + 2k 2 )k a (h)k a (l 2 ) + hkdikaih) + k a (l 2 )) 

m m Xa 
r r Xp + 7 
a = hk d (k a (l 2 ) - k a (h)) 

P = kMKik) + k a (h)) + (r + 2k 2 )k a (h)k a (l 2 ) 



We take the limit of these expressions as A — > oo: 
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We arrive at the new steady state: 

M = mk d^kik d + k a (l 2 )r 2 

1 r k x k d k s + k p r 2 
M 2 = ^ k Mr 2 + 2k d h 
r k d kik s + k p r 2 

A ^~( l + ttfft ) < R35 > 

r V kikdks + k p r 2 J 

m ( k\k d k \ 

^2 = -l- ,,, 1 " (B.36) 
r V k\kdk s + /c p r 2 / 

It is simple to verify that similarly to the original system where A ~ 0(1), the 

qualitative behavior is the same with respect to k\ and k 2 . If k± — 0, then the 

compartments reach the same steady state as when they were not connected. 

mk d k a (l 2 )r 2 mk d 



Mi = 



k p r 2 rk a (h) 
mk d 



M< 



2 



rk a (l 2 ) 

A X = A 2 = — 

r 

The case k 2 3> 1 and A 3> 1 is discussed in the main text. We conclude that regardless 
of the assumption we have make about A, we must have k\ ^> k 2 in order to have the 
desired dynamics. 

Finally, we examine the absolute difference between \A ± — A 2 \. As previously, we can 
assume without loss of generality that k 2 = 0. Then \Ai — A 2 \ is bounded below by 
zero and above by 

^ ^ 2m k\k d k 

t k\k d k s -\- vk p 

We consider A 1 — A 2 as a function of k and k s , as before. The same analysis as in 
the A 0(1) case shows that for a fixed concentration difference, k a (li) — k a (l 2 ), 
the optimal concentration range will be where kik d {k a {l\) + k a {l 2 )) = r. We have 
shown that for the particular cases we have considered, our equations have the same 
qualitative behavior for A ~ 0(1) and A> 1. This is sufficient for our purposes, but 
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we note that in order to make this statement rigorous, we would have to compare the 
approximate solution based on the separation of time scales (i.e. on the assumption 
that A is large) with the exact analytical solution. 

B.3 Analytical solution and approximation 



We consider the system of equations: 
dM 1 



dt 



dM 2 
~df 



m + \(-k a (h)M 1 + k d A x ) + h(M 2 - M ± ) 



dAi 
~dT 



-rA 1 + A(fc (Zi)M! - fc d Ai) 



m + \{-k a (l 2 )M 2 + k d A 2 ) - k 1 {M 2 - M x ) 



dA 2 
dt 



= -rA 2 + X(k a (l 2 )M 2 - k d A 2 ) 



MM = M 2 (0) = A 1 (0) = A 2 (0) = - 

Based on our previous analysis we assumed that the flux of A was much smaller than 
the flux of M, ki, so we set the rate of flux of A to be zero. We rewrite the our 
equations in matrix form. 



dy(t) 
dt 



with 



= Dy + h 

y(0) = m 



M x {t) 

M{t) 

M 2 (t) 

Mt) 
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D = 



Xk a (h) 

h 





\k d 
-{r + fed) 





h 
o 

-(Xk a (l 2 ) + h) 
Xk a (l 2 ) 






\k d 
-{r + k d ) 



h 



and 



2/o 



m 

m 


m fed 
r fe a (/ ) 

m 
r 

m _fcd_ 

r k a (l ) 



We want to find X, A such that -D = XAX 1 . Now our system becomes 

^ -^XAX^y + h (B.37) 
AX^ + X^/i (B.38) 



Define v — X l y and h = X 1 h, so 

^ A 7 

— = Aw + /i 
at 

v(0) = X-% 

By making one more substitution, and letting w — v + A^h we obtain 

dw 

— = Aw 
dt 

w{0) = v(0) + A~ l h 
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The solution to this equation is w = w(0)e At , and by making the appropriate substi- 
tutions again, this gives the solution to the equation IB. 381 to be 

y = (y + D- 1 h)e At -D- 1 h 

In order to find the explicit formula for y, we must find A, the diagonal matrix of 
eigenvalues of D and D~ x in terms of our parameters. In spite the fact that this is a 
problem with some symmetry, finding the eigenvalues and the inverse of the four-by- 
four matrix, D is difficult even with Matlab's Symbolic Math Toolbox. We leave the 
exact solution in this form. 

It is possible to approximate the exact solution to equations IB. 1 lllB. 151 in case A ^> 1. 
The solution is similar to the solution of the equations for perfect adaptation in 
Appendix IB. II 

Depending on the size of the flux k\ between the two compartments we can consider 
two cases. First, we assume that the flux between the two compartments is very fast, 
and in fact, k\ = 0(A). Based on our intuition developed by the steady state solution 
and the approximate solutions, we expect in this case the greatest change to be that 
the fast time scale, rj has to depend on both k a (li) and k a (l 2 ). If we write down the 
equations that apply on the fast time scale, 

^p- = -k a (h)M 1 + k d Ax + kx(M 2 - Mi) 
ar 

^ = k a (l 1 )M 1 - k d A, 
dr 

^ = -k a (l 2 )M 2 + k d A 2 - h(M 2 - MO 
ar 

^ = k a (l 2 )M 2 - k d A 2 

we see that now we must sum all four components to get a constant, i.e., Mi + 
M 2 + Ai + A 2 = C, so the four equations are coupled. In fact, solving this system of 
equations is not simpler than providing the exact solution, therefore we do not pursue 
this line of investigation. 
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Now we examine the case when k\ ~ 0(1), so on the fast time scale the it is still 
true that A\ + M\ = C\ for a constant C\, and similarly, A 2 + M 2 = C 2 . The same 
calculations we used in App endix IB . 1 1 apply, and we can obtain the solution to on the 
fast time scale: 

A x {r) = (A sl - A n )e- r ^ T + A n (B.39) 
M x {t) = (M i - M n )e- r ? lT + M n (B.40) 

A 2 (r) = (A s2 - A 12 )e- r ^ + A 12 (B.41) 
M 2 (r) = (M 02 - M 12 )e- r ^ T + M 12 (B.42) 

where we have defined the constants 

Tfi = k a (h) + k d , r f2 = k a (l 2 ) + k d 

mk d 



Moi = M 02 



r£; a (Z 



m A; d I + {k d /k a (l )) _ m k d 1 + (k d /k a (l )) 

11 ~ r fc a (Z x ) 1 + (k d /k a (h)) ' 12 " r fc a (/0 1 + (k d /k a (l 2 )) 

A A 7X1 

A sl = A s2 = — 
r 

_ ml + (k d /k a (l )) _ ml + (k d /k a (l )) 



On the slow time scale it remains true that A\ = 1 Mx , and A 2 = ^ M 2 . 

^d " ^d 



rl + (WW r l + (k d /k a (l 2 )) 

and A 2 = ^1. 

^d ^d 

As before, substituting these expressions into equations IB. 1 lllB~T5l we arrive at a new 
system of equations: 

dM 1 



m + k 1 (M 2 - Mi) (B.43) 
i^Mi (B.44) 



fc«(Zi) dMi fc a (Z 



-r- 



A; d (it fc d 
^ = m - jfei(Af 2 - Mi) (B.45) 
^i) dM = _ r ^) M2 

As before, we add equations IB. 431 and IB.44[ and equations IB.45I and IB.46I We arrive 
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at 



dM\ mk d 

dt k d + k a (h] 
dM 2 mkd 



+ 



+ 



k\k d 



k d + k a (h 
k\k d 



-{M 2 -M 1 )-r 



k a (h 



dt k d + k a (l 2 ) k d + k a (l 2 
The solution to eqns. IB.471 IB.48I is given by: 



(M 2 -M x )-r 



k d + k a (h) 

k a {l 2 ) 
k d + k a (l 2 ) 



Mi 

M 2 



where 



M(t) = (M(0) + D- l h)e M - D~ l h 



D 



-rk a (h)-kik d 
kd+k a (h) 

k\k d 
k d +k a (h) 



k\k d 



k d +k a (h) 
-rk a (h)—k-ik d 
k d +k a {l2) 



(B.47) 
(B.48) 



h 



mk d 



k d +k a (h) 

mk d 
k d +k a (h) 



M(0) 



c 2 



We define f3 = Tr(D) and 7 = det(-D). Then the eigenvalues Ai, X 2 of D can be 
found as follows: 

Ai,2 - 2 

By making the appropriate substitutions and carrying out the calculations, the dis- 



criminant \/ (5 2 — can be reduced to: 



(hkd + rkaih)) (hkd + rkaih))^ 



(hk d ) 2 



(k d + k a (h)) (k d + k a (l 2 )) J (k d + k a (h))(k d + k a (l 2 )) 

If we assume that the second term is much smaller than the first one, then the above 



expression greatly simplifies. This is true if 



(r + k^kaikaih) - k a (l 2 )) > k x k d 



(B.49) 
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Equation IB.49I implies that our approximation is appropriate when the ligand con- 
centrations in the two compartment are very different, i.e. when k a (l\) — k a (l 2 ) 3> 1- 
Now we are able to find the two eigenvalues: 

_ rk a (h) + hk d rk a (l 2 ) + k x k d 
1,2 ~ k d + k a (h) ' k d + k a (l 2 ) 

The two eigenvalues define the two slow time scales, 

rk a (h) + hkd 



si 



r S 2 



k d + k a (h) 
rk a (l 2 ) + hk d 



kd + k a (l 2 ) 

We notice that if the flux k\ is small, then we have recovered the slow time scale for 
each compartment independently of each other. This result confirms conclusions we 
have drawn from our steady state analysis. 
Next, we find D~ l h: 

m k d (rk a (l 2 ) + hk d ) 



d 1 = D~ l hi 



r rk a (li)k a (l 2 ) + hkdihih) + k a {l 2 )) 



a n-U m kdh(k d + k a (h)) 

a 2 — L) a 2 - 



r rk a (h)k a (l 2 ) + kik d (k a (h) + k a (l 2 )) 
The solution to eqns. is: 

Mi = (ci + rfi)e ( - rsl ' } - d l (B.50) 
M 2 = (c 2 + d 2 )e^ r ^ - d 2 (B.51) 

The equations for Mi and M 2 also determine the expressions for A\ and A 2 . 

A, = ^#(ci + dM~ Tslt) - k -^T^ (B.52) 
k d k d 

A 2 = M^(c 2 + d 2 )e^ - (B.53) 
k d k d 

Now we can match the fast and slow solutions as before to determine the constants 

c\ and c 2 , and to give arrive at the full approximation. Taking the limit as r — > 00 
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of the equations on the fast time scale, and setting this equal to the initial condition 
of the equations of the slow time scale gives us: 

d + d l - di = M n 
c 2 + d 2 - d 2 = M12 

— (Ci + di) : = An 

k d k d 

k a {h) , . , \ k a {l 2 )d 2 . 

— — (c 2 + d 2 ) : = A 12 

kd k d 

Thus the approximate solution to our system of equations 

Ai = (A sl - A n )e-^ xt + ^(A n + dje™ - (B.54) 
Mi = (Mm - Mn)e -r ' lA * + (Mn + di)e ( - rslt) - di (B.55) 



MM (Al2 + ( | 2)e M_W 2 

M 2 = (M 02 - Mi 2 )e-^ A * + (Mi 2 + d 2 )e^ r ^ - d 2 (B.57) 



(A s2 - A l2 )e-'^ + ^(A 12 + d 2 )e^ - (B.56) 

kd k d 



B.4 Calcium switch 

We return to the experimental observation that the cytosolic calcium concentration 
can change the turning behavior of a growth cone. Let us assume that a netrin-1 
gradient is set up outside the cell. Recall that in a cell with normal cytosolic calcium 
concentration a gradient develops with the high calcium concentrations facing the 
source of netrin-1, and the growth cone responds with attractive turning. However, in 
the same netrin-1 gradient a cell whose cytosolic calcium has been depleted before the 
trial responds with repulsive turning (eventhough the high calcium concentrations still 
face the source of netrin-1). Such a behavior is possible by making some assumptions 
regarding k a , the rate at which A is produced in our model. 

We assume that the production rate of A depends both on the ligand concentration 
and the cytosolic calcium concentration, k a (l,Ca). Further, we want k a to be such 
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that for large values of calcium ^ > 0, so k a is an increasing function of the ligand, 
and for small values of calcium ^ < 0. Without the constraints of experimental 
data, many such functions can be found. We chose k a = exp( n^jT^q^ ) where a, b 
and c are new constants, and Cab is the normal cytosolic calcium concentration. k a 
is illustrated in Figure IB. II 




Figure B.l: k a (l,Ca), the production rate of the modified substance, A. 1: ligand, netrin-1, Ca: 
cytosolic calcium concentration, a = 0.01, b = 1, c = 1, Cab = 0.2 

Using this rate we can simulate what happens with our system of equations in the 
same ligand gradient when the calcium level is above and when the calcium level is 
below the baseline, Cab- Changing the calcium level in the simulations corresponds 
to changing the value of the parameter Ca in the expression for k a (l). This implicitly 
implies that we take the calcium level to be spatially uniform inside the cell. Figure 
IB. 2 1 shows that in the normal cytosolic concentration we get a gradient of the adapted 
substance, A when a netrin-1 gradient is presented, with the higher concentration of A 
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corresponding to the higher concentration of I. (As before, the "gradient" of A is based 
on two values only, A 1 = A(l) in the left hand compartment and A 2 = A(2) in the 
right hand compartment.) Lowering the cytosolic concentration level, and presenting 
the cell the same ligand gradient results in a gradient of A where now the level of A 
is lower in the compartment corresponding to the higher ligand concentration. 



Two compartment moOel, k a (l,Ca| 




ca =0.5 ca=0.8 



Two compartment moOel, k a (l,Ca| 




ca b =0.5, ca=0.1 



1 1.1 1.2 1.3 1.4 



1.6 1.7 1.6 



1 1.1 1.2 1.3 1.4 



1.6 1.7 1.8 



Figure B.2: Numerical simulations of the two compartment model using k a — cxp( ''j^'^) ) ■ The 
first figure shows a cell where the cytosolic calcium level is above the baseline. The second figure 
shows a cell in the same ligand gradient where the calcium level is below the baseline. 



Both panels of the figure are run to the same time, t=50 seconds. The simulations 
are based on the same Matlab code as the previous two compartment models, except 
that now k a {l) is changed. 

This model presents a very simple explanation of how calcium levels influence turning 
behavior. Only the overall calcium concentration is considered by the simulations, 
the spatial and temporal gradients of calcium are not. Further investigations are 
necessary to create a more realistic description of the behavior. 
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Appendix C 

Sample Matlab code 

We are interested in investigating the deformations of two coupled Kelvin bodies in 
(i) steady flow, F = Fq and (ii) oscillatory flow, F = Focos(ut). Equations 14.91 and 
14.101 were solved with a Matlab code given below. The solution presented here is for 
oscillatory flow. 

First we can see the program twobodies_aF .m that asks for the input from the user 
and displays the results of the simulations. The user defines the coefficients for the two 
Kelvin bodies, and the program stores this information in the matrix B. In this code a 
function, parallel2_aF .m is called which actually computes the solution with a four- 
stage Runge-Kutta method. Then, as the output of parallel2_aF .m, the solution 
of the differential equation is returned to the matrix u_l in twobodies_aF .m, the 
solution (in this case u(t)) is plotted. 

r=input ( 'Two-body system with the bodies connected in parallel'); 
r=input (' Coefficients for ith body are: 
[\mu_{0i} \mu_{li} \eta_{li}] '); 

B=zeros(2,3) ; 

B(l ,: )= input ( 'Coefficients for Body 1:- '); 
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B(2, :)=input( 'Coefficients for Body 2:- '); 



h = 



0.1; 



% size of the time step 



N0=5; 



N = NO/h; 



x=[l:l:N] ; 



% length of a, u, F 



u_l=paralle!2_aF(B,h,N) ; 



°/„ function call 



ul=u_l(l,l:l:N) ; 
F=u_l(3,l:l:N) ; 
a=u_l(2,l:l:N) ; 
a=a./F; 

plot (x,u, 'r- . ' ) 

ylabel( 'u(t) , deformation') 

xlabel ( 'time , t') 

title('u vs t for different values of \eta_{12}, 
F=F_0 cos(\omega t)') 

Now we can look at the code for the actual ODE solver, parallel2_aF.m. Here 
only the main loop of the program is included which contains the fourth-order four- 
stage Runge-Kutta method. In the actual program the matrices Ml, M2, M3 and M are 
defined to be A, D, c and A~ l D, respectively. The full code contains the initialization 
of all the appropriate variables. This program also calls a function, ve2_aF.m, given 
below. The input of parallel2_aF .m function is the matrix of coefficients of the 
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Kelvin bodies, denoted by B , the step size, h, and the length of the solution vector, 
N. The output of this function is also a matrix, called ul_plot whose first row is u(t), 
second row is a(t)F(t), and third row is F(t). u(t) and a(t)F(t) are obtained as the 
solutions to the differential equations. F(t), the third row of the matrix, is simply 
F cos(ujt) evaluated for each time step. 

function f =parallel2_aF(B,h,N) ; 
for k=l:N 

F(k+l)=FO*cos(w*(t+h)) ; °/„ oscillatory 

kl=u; 

k2=u+ (1/2) *h*ve2_aF (M , inMl , M3 , kl , t+ (h/2) ) ; 
k3=u+ (1/2) *h*ve2_aF (M , inMl , M3 , k2 , t+ (h/2) ) ; 
k4=u+h*ve2_aF (M , inMl , M3 , k3 , t+h) ; 

u_new=u+h*((l/6)*ve2_aF(M,inMl,M3,kl,t) . . . 
+(l/3)*ve2_aF(M,inMl,M3,k2,t+(h/2)) . . . 
+(l/3)*ve2_aF(M,inMl,M3,k3,t+(h/2)) . . . 
+ ( 1/6) *ve2_aF (M , inMl , M3 , k4 , t+h) ) ; 

u = u_new; 
t = t+h; 

ul_plot(l:2,k+l)=u; 
ul_plot(3,k+l)=F(k+l) ; 



end 
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f =ul_plot ; 



Finally, we can look at ve2_aF.m. The input of this function is A~ X D, denoted here 
by M; A -1 , denoted by inMl; c, denoted by M3; a vector which consists of u(t) and 
a(t)F(t) at the previous time step, denoted by u, and finally, t, the current time step. 
The output of this function is the right hand side of equation 14.91 for the appropriate 
time step with the given type of flow. Here oscillatory flow is shown. 



function f = ve2_aF(M,inMl,M3,u,t) 

F0=1; 
w=2*pi ; 

F=F0*cos(w*t) ; 
dF=-FO*w*sin(w*t) ; 

C=zeros(2, 1) ; 

C = M3*dF; 

C(2,l) = F+C(2,l) ; 



f = M*u+inMl*C; 
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